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Abstract 

The  deflection  limitations  of  electrostatic  flexure-beam  actuators  are  well  known  [64]. 
Specifically,  as  the  beam  is  actuated  and  the  gap  traversed,  the  restoring  force  necessary 
for  equilibrium  increases  proportionally  with  the  displacement  to  first  order,  while  the 
electrostatic  actuating  force  increases  with  the  square  of  the  potential  difference  across 
the  gap,  as  well  as  the  inverse  square  of  the  gap.  Equilibrium,  and  thus  stable  open-loop 
voltage  control,  ceases  at  one-third  the  total  gap  distance,  leading  to  an  unstable  actuator 
snap-in.  A  Kalman  Filter  is  designed  with  an  appropriately  complex  state  dynamics 
model  to  estimate  actuator  deflection  accurately  given  voltage  input  and  capacitance 
measurements,  which  are  then  used  by  a  Linear  Quadratic  controller  to  generate  a  closed- 
loop  voltage  control  signal.  The  constraints  of  the  latter  are  designed  to  maximize  stable 
control  over  the  entire  gap.  The  design  and  simulation  of  the  Kalman  Filter  and 
controller  are  presented  and  discussed,  with  static  and  dynamic  responses  analyzed,  as 
applied  to  basic,  100  pm  by  100  pm  square,  flexure-beam-actuated  micromirrors 
fabricated  by  PolyMUMPs.  Successful  application  of  these  techniques  enables 
demonstration  of  smooth,  stable  deflections  of  50%  and  75%  of  the  gap. 
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LINEAR-QUADRATIC  CONTROL  OF  A  MEMS  MICROMIRROR  USING 

KALMAN  FILTERING 

1.  Introduction 


1.1.  Background 

In  the  past  decade,  the  field  of  MicroElectroMechanical  Systems  (MEMS)  has 
enjoyed  rapidly  accelerating  scientific  research  and  commercial  adoption  [63]. 
Advancements  in  fabrication  techniques  and  simulation  tools  have  enabled  greater  system 
refinement  and  understanding.  As  an  illustration  of  the  maturity  of  the  field,  MEMS- 
based  accelerometers,  such  as  Analog  Devices’s  ADXL-50,  have  been  universally 
adopted  since  2001  for  airbag  deployment  [63],  while  MEMS  pressure  sensors  are  now 
standard  features  in  automobiles  [63],  and  chips  containing  millions  of  MEMS 
micromirrors,  with  mean  times  between  failure  on  the  order  of  20  years  and  a  trillion 
cycles,  are  commercially  available  [1].  Single  micromirrors  with  two-dimensional 
scanning  capability,  as  those  shown  in  Figure  1,  are  commercially  available  “off-the- 


Figure  1:  Packaged  scanning  micromirrors, 
commercially  available  from  Adriatic  Research, 
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shelf’  in  standard  dual  inline  packages. 

Fabrication  variability  and  analytical  complexity,  however,  continue  to  limit  the 
growth  of  MEMS.  For  the  PolyMUMPs  process  in  particular  (albeit  a  pathological 
example  primarily  used  for  education),  parameter  deviation  of  20%  from  nominal  values 
published  by  MEMSCAP  has  been  reported  [5],  all  but  eliminating  the  possibility  of 
highly  accurate  modeling.  Exacerbating  this  problem,  variation  in  device  dimensions  can 
be  magnified  in  system-level  parameters  quadratically  or  even  cubically,  depending  on 
application.  Consider,  for  example,  a  rectangular  flexure  beam  of  the  classical  comb 
resonator;  a  20%  decrease  in  width  of  this  beam  will  result  in  a  50%  smaller  moment  of 
inertia,  which  will  be  directly  reflected  in  the  effective  spring  constant  of  the  beam,  and 
thus  the  resonant  frequency  of  the  device.  Combine  this  magnification  with  complex 
geometries  and  strong  coupling  between  the  applicable  physics  regimes,  and  traditional 
approaches  to  actuator  design  and  control  become  nearly  intractable. 

The  Kalman  Filter  (KF)  is  uniquely  suited  to  estimating  parameters  and  operating 
variables  that  are  otherwise  difficult  to  measure.  More  commonly  used  in  autopilot  and 
inertial  measurement  systems,  KFs  exploit  knowledge  about  system  and  sensor  dynamics, 
noise  sources,  and  initial  conditions  to  provide  a  running,  quantitative  characterization  of 
a  system.  In  this  capacity,  KFs  act  as  a  class  of  stochastic  observers,  which  can  then  be 
used  to  extend  the  capability  of  a  controller  by  increasing  the  observability  of  that  which 
is  to  be  controlled.  With  accurate  estimates  of  physical  parameters,  the  observer 
increases  the  controller’s  insight  into  the  system,  generally  allowing  more  efficient 
control  actions. 
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The  flexure-beam  micromirror  device  (FBMD),  shown  in  Figure  2,  exemplifies 
the  type  of  nonlinear,  highly  dynamic  system  for  which  a  KF  may  be  employed  to  great 
benefit.  A  well-known  artifact  of  most  electromechanical  microactuators  is  the  so-called 
pull-in  voltage;  that  is,  the  voltage  at  which  mechanical  restoring  forces  cease  to  balance 
Coulombic  attractive  forces,  and  the  electrodes  snap  together  as  a  result  of  the  lost 
equilibrium.  This  condition  uniformly  occurs  at  approximately  one-third  of  the  total 
available  travel  distance,  commonly  referred  to  as  the  “gap,”  leaving  two-thirds  without  a 
stable  operating  point.  Furthermore,  such  an  uncontrolled  impact  accelerates  contact 
wear  and  causes  hysteresis  in  the  gap  versus  applied  voltage,  since  stiction  (a  Van  der 
Waals  force  that  causes  two  surfaces  in  contact  to  be  “sticky”)  adds  to  the  force  needed  to 
return  to  equilibrium.  In  the  worst  case,  this  stiction  is  pennanent  and  leads  to  device 
failure,  which  constitutes  the  major  factor  in  determining  device  reliability. 

Only  for  testing  purposes,  however,  are  FBMDs  constructed  such  that  the  gap  is 
directly  measurable.  This  is  usually  accomplished  by  back-etching  through  the  bottom 
electrode  of  the  FBMD,  shining  a  laser  on  the  moving  plate,  and  collecting 
interferometric  data,  thereby  optically  measuring  the  distance  travelled  by  the 
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micromirror  [5].  Another  method  of  collecting  the  same  data  without  the  back-etch  is  a 
down-looking  interferometer,  such  as  the  popular  Zygo  white-light  interferometer; 
placing  a  control  sensor  in  the  path  of  incoming  signal,  however,  would  inevitably 
decrease  the  utility  of  the  micromirror.  Furthermore,  interferometry  takes  time  to  acquire 
data  and  signal  processing  to  convert  data  into  a  meaningful  measurement  of 
displacement,  which  severely  handicaps  the  bandwidth  of  any  controller  implementation. 
The  additional  fabrication  processing  steps,  additional  light  source,  and  associated 
electronics  needed  to  perfonn  this  measurement,  moreover,  add  system  complexity  so  as 
to  be  inappropriate  for  an  array  implementation,  as  well  as  cost  prohibitive  for 
implementation  on  any  commercial  scale.  The  work  presented  here  attempts  to 
demonstrate  the  utility  of  using  a  KF  to  estimate  this  gap  purely  as  a  function  of  input 
current  and  measured  voltage  across  the  FBMD,  using  known  system  dynamics,  initial 
conditions,  and  approximate  sensor  noise  strengths,  enabling  nearly  instantaneous 
measurement. 

1.2.  Related  Work 

Much  work  has  been  done  to  control  electrostatic  MEMS  actuators,  most  of 
which  with  the  goals  of  extending  the  operating  travel  range,  and/or  increasing 
positioning  accuracy  [5],  These  efforts  can  be  divided  into  several  categories:  1) 
geometrical,  in  which  the  device  structure  is  adjusted  to  increase  the  effective  mechanical 
spring  constant,  or  increase  the  stroke  length;  2)  open-loop  control,  which  manipulates 
the  input  signal  to  create  desired  effects  based  on  empirically  or  analytically  identified 
system  parameters;  and  3)  closed-loop  control,  which  uses  some  sort  of  error  signal  to 
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tune  the  input  signal.  With  respect  to  an  intended  application,  each  method  has 
advantages  and  disadvantages. 

Changing  the  geometry  of  the  device  often  complicates  the  fabrication  process, 
with  the  benefit  of  relatively  simple  operation.  This  is  usually  done  to  enhance  the 
effective  spring  constant,  or  create  and/or  control  a  nonlinear  restoring  force,  thereby 
extending  the  actuator  stroke.  Common  themes  include  a  phased  flexure  [6] [33],  such 
that  supports  and/or  beams  are  activated  as  deflection  increases,  effecting  a  nonlinearly 
increasing  spring  constant  with  increasing  deflection;  bottom  electrode  sizing  and 
positioning  [1 1][14],  in  which  the  capacitance  area  or  fulcrum  position  is  varied  such  that 
increasing  deflection  exposes  more  (or  less)  of  the  bottom  electrode;  and  optimization  of 
the  device  structure  with  respect  to  the  pull-in  voltage  [18]. 

Open-loop  control  attempts  to  maintain  an  unstable  position  by  varying  the  input 
waveform,  despite  not  knowing  the  real-time  actuator  deflection.  This  variation  is  based 
upon  insight  into  the  system  dynamics,  usually  from  analytical  modeling,  or  empirical 
system  identification.  Waveforms  may  take  the  form  of  pulsed  voltage  [13] [29]  [35],  or 
charge  [7] [28],  or  of  continuous,  “preshaped”  signals  [12].  While  the  structure  itself 
remains  relatively  simple  in  terms  of  fabrication,  open-loop  control  requires  drive 
electronics  complex  enough  to  generate  non-trivial  wavefonns,  and,  more  importantly, 
generally  features  less  accurate  deflections  than  similar  closed-loop  methods  as  a  result  of 
a  lack  of  robustness  to  process  variations,  environmental  challenges,  or  deviations  from 
system  concept  of  operations.  Inaccurate  deflection  control  does  not  prevent  snap-in,  but 
rather  modestly  extend  the  deflection  range. 
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Comprising  the  majority  of  research  in  the  area  is  that  of  closed-loop  control. 
Efforts  range  from  applying  classical  linear  controllers  to  simplified  equations  of  motion 
[3] [4]  [5]  [  1 0]  [  1 6]  [22]  [25]  [27]  [3 1  ]  [32]  [34]  [43] ,  to  adding  compensators  to  effect  feedback 
linearization  [8][17][2 1] [22] [28] [37] [38] [39] [40] [4 1  ],  to  generating  fully  nonlinear 
control  laws  [20]  [45]  [46]  [47]  [48][49].  Some  designs  employ  full-  or  reduced-order 
observers  [2]  [9]  [  1 5]  [  1 9]  [23]  [25]  [26]  [30]  [32]  [42]  [46],  which  estimate,  in  real-time,  the 
position  and/or  velocity  of  the  actuator  based  upon  measured  parameters  (e.g., 
capacitance).  These  observers  then  update  the  assumptions  made  by  the  controller, 
making  it  adaptive.  Less  common  methods  include  neural  networks  [1],  fuzzy  logic 
[18][19],  sliding  mode  control  [36],  port-controlled  Hamiltonian  systems  [22],  and 
passivity-based  control  [23]  [24] [26]. 

1.3.  Problem  Statement 

The  instability  of  electrostatic  flexure-beam  actuators  beyond  one-third  of  the  gap 
across  which  the  potential  is  applied  leads  to  a  nonlinear  “snap-in”  effect  that  limits  the 
effective  range  of  controllable  actuation  and  dramatically  reduces  operational  lifetime. 
From  a  control  design  standpoint,  this  problem  is  exacerbated  by  limited  system 
observability  (nominal  structures  exclude  feedback  sensors)  and  wide  parameter 
variability.  The  Kalman  Filter,  used  in  conjunction  with  a  Linear-Quadratic-Gaussian 
(LQG)  controller,  is  uniquely  suited  to  estimate  unobservable  states  in  the  presence  of 
such  parameter  variability  and  noise  sources,  owing  to  its  simple  measurement  system, 
scalability  to  large  arrays,  and  straightforward  digital  implementation.  This  powerful 
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combination  enables  real-world,  full-gap  actuator  positioning  and  controllable  snap-in, 
thereby  eliminating  the  dominant  failure  mechanism. 

1.4.  Scope 

This  work  is  limited  to  a  proof-of-concept  observer/controller  pair  as  applied  to  a 
single  micromirror  geometry.  While  the  approach  presented  is  equally  valid  for 
alternative  geometries,  material  systems,  and  actuator  schemes,  only  one  will  be 
explored.  Furthermore,  no  attempt  to  create  a  deployable  system  will  be  made.  Issues  to 
be  considered  for  robustness  enhancement,  such  as  external  shock  or  vibration  forces, 
changes  in  ambient  temperature  or  pressure,  and  performance  degradation  over  time,  will 
be  left  to  follow-on  research.  Lastly,  since  this  work  is  a  proof-of-concept,  hardware 
implementation  will  be  simulated  and  limited  to  that  which  is  required  for  developmental 
validation  of  the  control  algorithms;  electronics  integration,  footprint  minimization,  and 
packaging  will  also  be  left  for  future  work. 

1.5.  Preview 

This  research  is  presented  in  five  chapters.  Chapter  1  introduces  the  problem, 
motivates  the  solution,  and  delineates  past  research  with  similar  aims.  The  problem 
statement  clearly  identifies  the  goal  of  the  work,  while  the  scope  sets  boundaries  and 
specifies  starting  assumptions. 

Chapter  2  presents  background  theory  in  sufficient  detail  such  that  the  reader  can 
understand  the  design  process  described  in  later  chapters.  In  particular,  a  summary  is 
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included  of  PolyMUMPs  fabrication  steps,  the  mechanical  and  electrostatic  physics  of 
flexure-beam  actuators,  and  Kalman  Filter  and  LQG  controller  operation. 

Chapter  3  describes  the  modeling  and  the  control  law  development.  The 
COMSOL  Multiphysics  modeling  and  simulation  software  is  used  to  generate  nominal 
charge  and  position  versus  time  trajectories,  upon  which  the  KF  is  based.  Lastly,  the 
LQG  controller  is  derived  from  actuator  position  and  velocity  requirements. 

Chapter  4  discusses  the  perfonnance  of  the  control  laws  developed  in  Chapter  3. 
Results  from  MATLAB  simulations  are  compared  to  the  nominal  trajectories  by 
statistical  analysis,  with  the  system  perfonnance  indicated  by  the  state  error  means  and 
covariances.  A  microminor  fabricated  in  PolyMUMPs  is  then  used  to  accomplish  a 
hardware-in-the-loop  test  for  real-world  perfonnance,  the  latter  characterized  by 
comparing  discrete  gap  measurements  to  the  KF  estimate. 

Chapter  5  provides  conclusions  and  suggestions  for  future  work.  The  system  is 
judged  based  upon  the  results  in  Chapter  4,  and  perfonnance  shortfalls  highlighted  and 
explained. 
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2.  Background 


2.1.  Chapter  Overview 

The  purpose  of  this  chapter  is  to  summarize  the  foundations  upon  which  this  work 
is  built.  First,  the  PolyMUMPs  fabrication  process  is  described  step  by  step,  and  key 
artifacts  highlighted.  Second,  the  FBMD  layout  is  presented  in  the  context  of 
PolyMUMPs  fabrication.  Third,  a  first-order  analytical  model  is  derived  from  the 
respective  physical  regimes  (mechanical,  electrostatic,  and  fluid  dynamics).  After  this 
derivation,  the  stability  of  the  model  is  analyzed  for  controller  suitability.  Last,  the 
Kalman  Filter  and  Linear  Quadratic  Gaussian  control  are  introduced  for  further 
discussion  in  Chapter  4. 

2.2.  Fabrication 

The  MEMSCAP  Multi-User  MEMS  Processes  (MUMPs)  is  employed  to  fabricate 
prototype  MEMS  devices  for  government,  industry,  and  academia  worldwide.  In 
particular,  PolyMUMPs  features  three,  surface-micromachined  polysilicon  layers,  two  of 
which  are  releasable  (that  is,  the  layer  immediately  beneath  can  be  etched  away  to 
“release”  the  polysilicon  above  it),  and  one  metallization  layer.  The  fabrication  process 
is  fixed  and  enforced  by  design  rules  specified  by  MEMSCAP  [66].  These  are  described 
below  and  illustrated  in  Figure  3  in  order  to  elucidate  eccentricities  of  the  process. 

First,  n-type  (100)  silicon  wafers,  150  millimeters  in  diameter,  are  pre-treated  by 
heavily  doping  the  surface  with  phosphorus  in  order  to  help  prevent  charge  accumulation 
between  the  substrate  and  isolation  layers.  This  isolation  layer  is  then  created  by  0.6  pm 
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Figure  3:  PolyMUMPs  process  layers  [66] 


of  low-stress  low-pressure  chemical  vapor  deposition  (LPCVD)  silicon  nitride,  followed 
by  a  0.5  pm  polysilicon  layer,  called  polyO.  This  polysilicon  layer  is  also  deposited  by 
LPCVD,  photolithographically  patterned,  and  etched  via  plasma  etching.  The  first 
sacrificial  layer,  oxide  1,  consisting  of  2.0  pm  of  phosphosilicate  glass  (PSG)  [66],  is 
deposited  next  by  LPCVD  at  580  °C  and  then  annealed  at  1050  °C  in  argon.  This  anneal 
reduces  residual  stress  and  dopes  the  lower  polysilicon  layer  to  a  concentration  of  lxl 019 
cm'  [67].  Stiction-preventing  dimples  are  created  by  reactive-ion-etching  (RIE)  0.75  pm 
holes  into  the  PSG  layer  at  this  stage.  Anchor  points,  which  connect  the  second 
polysilicon  layer  and  the  substrate  {anchor  1),  are  similarly  produced  by  RIE  immediately 
following  the  dimple  etch. 

Starting  with  the  deposition  of  the  polysilicon  layer,  the  process  is  essentially 
repeated,  with  some  differences.  The  second  polysilicon  layer,  polyl,  is  2.0  pm  thick  and 
deposited  by  LPCVD  at  this  step.  This  is  followed  by  a  PSG  layer  to  act  as  a  hard  mask 
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for  the  polysilicon  etch,  which  is  subsequently  removed  by  RIE.  Next,  yet  another  PSG 
layer  is  deposited  0.75  pm  thick  and  annealed  to  act  as  the  second  sacrificial  oxide  layer, 
oxide2.  Two  different  etch  masks  are  used  at  this  point:  the  first  to  create  a  connection 
between  the  second  and  third  polysilicon  layers  (polyl-poly2  via);  and  the  second  to 
connect  either  the  first  and  third  polysilicon  layers,  or  poIy2  to  the  nitride  layer 
( anchor2 ). 

The  last  polysilicon  layer,  po/y2,  is  now  deposited  1.5  pm  thick.  A  0.2  pm  PSG 
layer  is  again  used  as  a  hard  mask  and  dopant  source  for  the  third  polysilicon  layer;  each 
are  patterned,  and  the  PSG  removed  as  before.  The  final  layer,  metal  (0.5  pm  of  gold 
with  a  thin  chromium  adhesion  layer),  is  lithographically  patterned,  deposited  by  electron 
beam  evaporation,  and  patterned  using  lift-off.  With  the  fabrication  complete,  the  wafers 
are  diced,  sorted,  and  shipped. 

Three  points  must  be  emphasized.  First,  the  process  inherently  produces 
confonnal  layers;  that  is  to  say,  if  a  design  contains  a  feature  in  polyO,  the  poly  layers 
above  it  will  drape  over  the  polyO  shape,  like  a  rug  lying  over  a  book.  The  effect  is  that 
the  flexure -beam  actuators  fabricated  in  PolyMUMPs  are  not  the  straight  beams 
commonly  modeled  by  Newton-Euler  analytical  beam-bending  equations.  They  instead 
have  angles,  which  will  affect  the  restoring  force  generated  by  the  beam  for  a  given  load. 
Second,  successive  layers  of  polysilicon  undergo  fewer  anneal  stages,  which  affects  both 
the  conductivity  and  the  strength  of  the  layer.  In  particular,  the  Young’s  Modulus  of  the 
first  polysilicon  layer  has  been  measured  to  be  approximately  20%  less  than  the  second, 
while  exhibiting  half  the  resistivity  [68],  Third,  residual  stress  must  be  managed 
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thoughtfully  in  PolyMUMPS  devices.  As  an  artifact  of  bonding  two  materials  with 
different  coefficients  of  thermal  expansion,  residual  stress  manifests  most  apparently  as 
bowing  in  layers  of  poIy2  and  metal.  A  bowed  flexure  beam  will  deflect  differently  than 
an  unbowed  beam  of  similar  composition,  while  a  bowed  mirror  will  feature  a  different 
capacitance  due  to  the  non-uniform  gap.  These  effects  will  be  studied  in  Chapter  3. 


2.3.  Flexure-Beam  Micromirror 

The  micromirror  design  on  which  the  present  work  is  composed  is  shown  in 
Figure  4,  as  fabricated  in  PolyMUMPs,  described  above.  The  mirror  plate  is  a  100  pm 
by  100  pm  polyl  square,  with  five  etch  holes  to  ensure  all  oxide  is  removed  with  the 


polyl  Contact 


polyO  Contact 


Figure  4:  FBMD  layout 
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release  process  (shown  as  orange  in  Figure  4),  each  3  pm  by  3  pm.  A  polyO  square  of 
dimensions  equal  to  the  polyl  mirror  plate  is  immediately  below  the  plate  (orange  in  the 
figure)  and  serves  as  the  bottom  electrode,  typically  held  at  ground.  Note  that  the  polyO 
connection  between  the  bonding  pad  and  the  electrode  is  reflected  on  each  side  of  the 
square;  this  is  to  ensure  symmetry  in  the  three-dimensional  shape  of  the  flexures  given 
the  conformal  process.  That  is  to  say,  the  released  flexure  beam  will  have  upward  and 
downward  bends  with  respect  to  the  plane  of  the  substrate,  as  it  conforms  over  the  polyO 
connection.  The  symmetry  of  the  electrode  encourages  uniform  response  from  each  of 
the  four  flexure  beams.  The  flexures  are  each  100  pm  long  by  3  pm  wide  polyl  beams, 
attached  to  the  mirror  plate  by  3  pm  long  by  8  pm  wide  polyl  connectors.  These 
connectors  in  this  application  act  as  torsion  springs,  as  the  faces  attached  to  the  mirror 
and  flexure  beam  rotate  with  respect  to  each  other  as  the  mirror  deflects.  Located  on  the 
substrate  side  of  each  of  these  connectors  is  a  dimple,  which  prevents  the  mirror  plate 
from  coming  into  physical  contact  with  the  bottom  electrode,  minimizing  the  possibility 
of  stiction,  and  electrical  shorting.  10  pm  by  10  pm  polyl  anchors  attach  each  of  these 
beams  to  the  substrate,  thereby  fixing  one  end  to  a  rigid  post  (the  other  end  is  only  fixed 
to  the  torsion  spring  connector,  and  is  allowed  to  move  in  space).  One  anchor  point  is 
attached  to  a  50  pm  by  50  pm  poly 2  bond  pad,  which  creates  a  conductive  path  between 
the  mirror  plate  and  the  bond  pad  and  enables  an  external  voltage  source  to  generate  a 
potential  difference  between  the  plate  and  bottom  electrode. 
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Figure  5:  FBMD  modeled  mechanically  as  a  simple  harmonic  oscillator  and 
electrostatically  as  a  parallel  plate  capacitor 

2.4.  Analytical  Model 

Analytically  modeling  the  FBMD  promotes  understanding  of  the  pull-in 
phenomenon  and  allows  meaningful  exploration  of  possible  solutions.  As  mentioned 
above,  the  FBMD  is  held  in  equilibrium  by  two  opposing  forces:  the  mirror  plate  is 
electrostatically  attracted  down  towards  the  bottom  electrode  as  a  result  of  an  applied 
potential;  and  the  flexure  beams  mechanically  resist  the  electrostatic  force  and  restore  the 
plate  to  an  initial,  quiescent  distance,  go,  away  from  the  bottom  electrode,  as  illustrated  in 
Figure  5.  The  output  restoring  force  of  the  flexures  is  generally  modeled,  to  first  order, 
by  Hooke’s  Law,  i.e.,  linearly  proportional  to  the  amount  of  induced  deflection  in  the 
beams,  with  the  proportionality  constant  thought  of  and  referred  to  as  the  spring  constant, 
k: 
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Fm  =  kd  =  Hso  -  g) 


(1) 


where  go  is  the  quiescent  gap  distance;  d  is  the  mirror  plate  deflection  distance  from  go ; 
and  g  is  the  gap  distance.  Various  schemes  are  used  to  adjust  the  effective  mechanical 
restoring  force  in  order  to  meet  system  design  requirements  for  resonance  frequency, 
operating  voltage,  and  the  gap-versus-voltage  hysteresis.  These  schemes  include 
changing  the  beam  material  or  geometry  (e.g.,  I-beam-like  cross  section);  increasing  the 
number  of  flexures  or  attachment  point  (e.g.,  attaching  at  the  middle  of  the  mirror  edge 
rather  than  the  corner  shortens  the  beam  and  increases  k,  as  shown  in  Section  2.5.3);  and 
lever  ann  implementation  [14].  A  Duffing  spring  model  may  increase  accuracy  by 
adding  a  cubic  deflection  term  fitted  to  data  by  proportionality  constant.  This  model  is 
one  method  to  account  for  nonlinear  spring  softening  with  increasing  deflection.  The 
literature  reports  yet  more  complicated  effects  [11]. 

By  contrast,  the  electrostatic  force  is  inversely  proportional  to  the  square  of  the 
distance  between  the  mirror  plate  and  the  bottom  electrode  plate,  g,  and  proportional  to 
the  square  of  the  potential  difference  across  the  plates.  This  force  is  derived  from  the 
well-known  expression  for  the  energy  U  stored  in  a  parallel-plate  capacitor  in  steady-state 
and  the  definition  of  capacitance: 


U  = 


I^dr2 

2  g 


(2) 


where  Q  is  the  total  charge  on  the  two  plates;  C  is  the  capacitance;  V  is  the  potential 
difference  across  the  plates;  g  is  the  distance  between  the  mirror  plate  and  the  bottom 
electrode  as  above;  A  is  the  area  of  the  plate;  and  s0  is  the  free  space  permittivity.  The 
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force  is  then  defined  as  the  negative  gradient  of  the  potential  energy.  Since  this  system 
has  only  one  degree  of  freedom,  the  gradient  is  simply  the  first  partial  derivative  with 
respect  to  g: 


^  =  IMr  2 

%  2  g2 


(3) 


Although  Equation  (3)  increases  the  polynomial  order  of  the  static  equation,  as  well  as 
numerical  complexity  as  a  result  of  the  discontinuity  at  zero  gap,  it  is  far  from  a  complete 
characterization  of  reality.  First,  each  of  the  four  flexures  accumulates  charge,  thereby 
acting  as  four  additional  capacitors  with  respect  to  the  ground  plate  and  varying  with 
deflection.  Second,  the  parallel  plate  model  assumes  parallel,  infinite  plates.  Since  the 
FBMD  is  finite,  the  electric  field  is  not  wholly  contained  between  the  plates,  but  rather 
extends  outside,  giving  rise  to  what  are  known  as  fringing  fields.  Moreover,  the  parallel 
plate  model  neglects  the  thickness  of  the  plate;  in  reality,  the  sides  of  the  plate  also 
generate  electric  fields,  which  must  be  accounted  in  the  aggregate  electrical  potential 
energy.  Last,  the  parallel  plate  model  assumes  a  uniform  gap,  but  as  mentioned  in 
Section  2.2,  the  mirror  plate  is  not  perfectly  rigid  and  may  demonstrate  some  measure  of 
bowing  in  the  center.  Although  to  lesser  effect,  the  flexure  beams  may  also  demonstrate 
bowing  in  addition  to  the  bump  discussed  in  Section  2.2  and  shown  in  Figure  8,  either 
down  as  a  result  of  the  weight  of  the  mirror  plate,  or  up  as  a  result  of  residual  stress.  The 
capacitance  of  the  FBMD  is  simulated  in  Chapter  3,  and  the  magnitude  of  these  non-ideal 
extensions  will  be  quantified  to  show  that  the  fringing  fields  are  by  far  the  largest  non¬ 
ideal  effect;  flexure  capacitance  and  plate  thickness  contributions  are  negligible. 
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2.5.  Beam  Theory 

2.5.1.  Introduction 

Beam  deflection  is  described  analytically  by  the  Euler-Bernoulli  equation  for 
beam  bending.  As  a  differential  equation,  any  of  unifonn  or  point  loads  and  moments 
can  be  used  as  forcing  functions,  and  fixed  (clamped),  free,  or  simply  supported  end  faces 
may  be  used  as  boundary  conditions.  The  equation  is  a  simplification  of  linear  elasticity 
theory  with  the  following  assumptions:  the  beam  is  subject  only  to  pure  bending,  i.e.,  no 
torsional  or  axial  loads;  the  material  is  isotropic  and  homogeneous,  i.e.,  the  flexural 
rigidity  is  constant;  the  material  is  linearly  elastic  and  will  not  reach  the  plastic 
deformation  limit,  i.e.,  Hooke's  Law  is  obeyed;  the  beam  is  initially  straight  with  constant 
cross-section  throughout;  an  axis  of  symmetry  is  in  the  plane  of  bending;  the  proportions 
of  the  beam  are  such  that  it  would  fail  first  and  foremost  by  bending;  and  cross-sections 
of  the  beam  remain  planar  during  bending.  Ineluctable  deviations  from  these 
assumptions  are  addressed  after  the  derivation. 

2.5.2.  Derivation 

The  static  Euler-Bernoulli  beam  equation  may  be  shown  to  be  the  result  of 
combining  four  basic  relationships.  First,  a  kinematics  equation  specifies  how  the  beam 
moves.  In  one-dimensional,  linear  beam  theory,  this  amounts  to  describing  how  each 
point  in  a  lengthwise  cross-section  of  the  beam  is  displaced  with  deflection;  this 
displacement  is  equivalent  to  the  strain  in  the  beam.  By  assuming  that  deflections  are 
small  and  that  the  neutral  plane  does  not  change  in  length  under  load,  the  beam  bends 
into  an  arc  of  curvature  x,  and  the  angle  6  through  which  the  widthwise  cross-section 
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moves  can  be  equated  using  the  small  angle  approximation  to  the  negative  of  the  slope  of 
the  deflection  w: 


—  (4) 

Second,  a  constitutive  equation  relating  how  the  beam  moves  in  response  to  external 
forces  is  specified.  For  approximately  linear  materials,  Hooke's  Law  is  employed: 

(5) 

where  a  is  stress,  s  is  strain,  and  E  is  Young’s  Modulus.  Third,  using  so-called  force 
resultant  equations,  the  point-by-point  aggregate  effect  of  these  external  forces  in  a 
widthwise  cross-section  is  quantified  by  integrating  the  appropriate  stress  over  the  cross- 
section  and  equating  the  result  to  moments  M  and  shear  forces  V: 

(6) 


Last,  equilibrium  is  established  for  each  infinitesimal  length  by  equating  the  change  in 
shear  force  to  pressure  load  p  and  the  change  in  moment  to  the  shear  force  resultant: 

—  (V) 

—  (8) 

Using  algebra,  the  ratio  of  pressure  load  (force  per  unit  length)  to  rigidity  can  be  shown 
to  equal  the  fourth  derivative  of  deflection  with  respect  to  length,  i.e.,  the  canonical 
Euler-Bemoulli  beam  bending  equation: 
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(9) 


Each  successive  derivative  of  the  deflection  given  by  the  Euler-Bernoulli  equation 
has  a  corresponding  physical  interpretation.  The  first  derivative  with  respect  to  the  length 
along  which  the  beam  deflects  is,  for  small  values,  the  angle  between  the  neutral  axis  and 
the  beam.  The  second  derivative  is  the  net  moment  on  the  beam,  while  the  third 
derivative  is  the  net  shear  force.  Net  load  per  unit  length  is  represented  by  the  fourth 
derivative.  Each  of  these  may  balance  a  static  or  dynamic  forcing  term  that  models  the 
type  of  operation  the  beam  is  performing;  e.g.,  a  diving  board  would  have  a  point  load 
(force  multiplied  by  a  Dirac  delta  function  with  domain  shift  corresponding  to  the 
position  of  the  force)  as  a  forcing  term  for  the  fourth  derivative.  As  stated  above,  several 
support  types  can  be  represented  through  the  appropriate  use  of  boundary  conditions.  For 
example,  setting  the  deflection  and  slope  of  one  end  to  zero  models  a  fixed  support  at  that 
end,  while  no  deflection  and  no  net  moment  represents  a  pin  connection.  Boundary 
conditions  must  be  set  carefully  in  order  to  properly  capture  the  physics. 

2.5.3.  Lumped-model  Effective  Spring  Constant 

A  lumped-model  parameter  consolidates  the  gamut  of  complex  physical  processes 
into  a  “black  box,”  a  computationally  simple — or  at  least  more  straightforward — 
abstraction  that  approximates  an  output  for  a  given  input.  This  abstraction  often  imitates 
the  functional  fonn  of  a  more  familiar  relationship,  e.g.,  a  mass-spring  system  or  basic 
circuit.  In  the  present  case,  the  deflection  of  a  flexure  beam  in  an  FBMD  is  modeled  as  a 
simple,  linear  spring  obeying  Hooke’s  Law,  with  a  restorative  force  linearly  proportional 
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to  the  distance  deflected.  The  process  begins  with  the  Euler-Bemoulli  equation  for  beam 
supporting  a  point  load  at  the  free  end,  with  the  latter  realized  in  boundary  conditions.  In 
particular,  the  fourth  derivative  of  deflection  with  respect  to  beam  length,  i.e.,  pressure 
per  unit  length,  is  set  to  zero.  Assuming  a  constant  flexural  rigidity  El,  both  sides  of  this 
equation  are  integrated  four  times,  yielding  a  cubic  polynomial  for  deflection  versus 
length  position  with  four  indeterminate  constants  of  integration.  Boundary  conditions  are 
then  applied:  position  and  slope  of  the  fixed  end  of  the  beam  are  set  to  zero;  the  slope  of 
the  free  end  is  set  to  zero  as  a  result  of  being  attached  through  a  hinge  to  the  mirror  plate; 
and  the  shear  force  at  the  free  end  of  the  beam  is  set  to  equal  the  point  load.  After 
solving  for  the  deflection  at  the  free  end,  the  point  load  is  solved  for  as  the  dependent 
variable  in  terms  of  this  deflection;  the  result  is  an  equation  that  emulates  Hooke’s  Law 
for  linear  springs: 

- - (10) 

The  proportionality  constant,  twelve  times  ratio  of  the  beam  flexural  rigidity  and  the 
length  cubed,  is  considered  as  the  effective  spring  constant,  “lumped”  into  which  are  the 
essential  physics  quantified  by  the  Euler-Bemoulli  equation.  As  each  of  the  four  flexures 
contributes  this  restoring  force,  the  total  mechanical  restoring  force  is  provided  by  four 
springs,  thus  the  effective  spring  constant  is  multiplied  by  four. 

2.5.4.  Non-idealities  /  Extensions 

The  above  litany  of  constraining  assumptions  can,  in  several  cases,  be  modified  to 
extend  the  validity  of  the  equation.  First,  beam  dynamics  may  be  analyzed  by  the 
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addition  of  the  second  derivative  of  deflection  with  respect  to  time,  scaled  by  mass  per 
unit  length  of  the  beam,  to  the  left-hand  side  of  Equation  (9).  Second,  superposition  may 
be  employed  to  model  three-dimensional,  distributed,  transverse  loading,  as  well  as 
composite  beams  through  the  flexural  rigidity  term.  Materials  not  obeying  Hooke's  Law 
may  be  modeled  by  the  appropriate  constitutive  equation  describing  the  relationship 
between  stress  and  strain  of  the  material  system.  In  this  way,  viscoelastic  or  plastic 
deformations  and  nonlinear  material  behavior  are  incorporated,  thereby  generalizing  the 
Euler-Bemoulli  equation.  Geometrically  nonlinear  beams  may  be  accounted  for  by 
dividing  the  second  derivative  of  deflection  by  the  three-halves  power  of  the  slope 
squared  with  unity  offset.  This  divisor  enables  a  more  accurate  description  of  an  initially 
curved  beam,  e.g.,  a  cantilever  with  residual  stress,  as  physically  motivated  in  Section 
2.2.  Large  deflections  (i.e.,  bending  radius  equal  to  or  smaller  than  one-tenth  of  the 
cross-section)  may  be  approximated  by  multiplying  the  moment  of  inertia  by  a  function 
that  increases  inversely  with  the  radius  of  curvature,  and  by  adding  another  moment  term 
to  the  second  derivative  of  deflection  that  does  the  same.  Last,  thick  beams,  for  which 
the  transverse  shear  strain  is  non-negligible,  must  depart  from  Euler-Bemoulli  treatment 
entirely  and  be  analyzed  with  the  Timoshenko  beam  theory. 

2.6.  Squeeze-Film  Damping 

To  guarantee  a  thorough  analysis,  the  application  of  MEMS  devices  outside  a 
vacuum  should  take  into  account  the  effects  of  submersion  in  a  viscous  fluid.  For  the 
dynamic  operation  of  the  FBMD  used  in  the  present  work,  air  between  the  bottom 
electrode  and  mirror  plate  is  forced  outward  as  the  mirror  is  pulled  down,  while  it  is 
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conversely  pulled  inward  as  the  mirror  is  restored  to  its  initial  position.  This  movement 
of  air  dissipates  energy  in  the  system  by  opposing  the  motion  of  the  mirror  in  both 
directions.  As  such,  it  can  be  modeled  as  a  dashpot  on  its  effect  in  a  mass-spring  system. 

A  number  of  simplifying  assumptions  can  make  tractable  the  Navier-Stokes 
equation — a  nonlinear  partial  differential  equation  governing  fluid  dynamics.  First,  the 
fluid  is  assumed  to  be  isothermal  and  Newtonian;  that  is,  the  ratio  of  the  shear  stress 
exerted  by  the  fluid  (drag)  to  the  rate  of  strain  is  equal  to  the  viscosity  of  the  fluid,  a 
constant.  Ignoring  thermal  variation  is  equivalent  to  assuming  unifonn  proportionality 
between  density  and  pressure.  Second,  viscous  forces  (the  fluid’s  resistance  to 
defonnation  by  shear  or  tensile  stress)  are  assumed  to  dominate  over  inertial  forces  (the 
fluid’s  resistance  to  changes  in  momentum),  known  as  Stokes  flow,  or  creeping  flow;  this 
is  due  to  the  dimensions  of  the  FBMD  being  small  enough  that  the  dynamic  viscosity  of 
air  is  much  larger  than  the  mass  of  the  air  in  the  cavity.  As  a  result,  the  general  Navier- 
Stokes  equation  can  be  simplified  to  the  Reynolds  equation,  another  second-order  partial 
differential  equation  closely  resembling  the  classical  heat  equation  with  internal  heat 
generation  [63].  Further,  assuming  a  unifonn  fluid  thickness,  small  pressure  variation 
with  respect  to  the  ambient,  and  small  mirror  plate  displacements  reduces  the  Reynolds 
equation  to  a  simple  Poisson’s  equation,  readily  solved  analytically  using  a  Green’s 
function.  The  result  is  in  a  form  germane  to  the  present  work,  namely,  a  function  that  is 
linearly  proportional  to  the  velocity  of  the  mirror  plate  to  within  a  constant. 

A  more  exact  formulation  takes  into  account  unique  features  of  the  MEMS,  such 
as  compressibility  effects,  slip-conditions,  and  large  mirror  plate  displacements. 


22 


Compressibility  effects  describe  the  disproportionate  fluid  outflow  as  a  result  of  mirror 
displacement  and  can  be  treated  as  additional  stiffness  in  the  mass-spring  system  [51]. 
Slip-conditions  reduce  the  damping  constant  by  taking  into  account  the  lack  of  continuum 
as  the  mean  free  path  of  fluid  molecules  becomes  significant  with  respect  to  the  fluid 
thickness  [56].  The  solution  to  the  Poisson’s  equation  displays  a  strong  dependence  on 
film  thickness,  so  it  is  intuitive  to  understand  that,  as  the  mirror  plate  displacement 
increases  with  respect  to  the  nominal,  the  damping  coefficient  will  increase  as  well  [69]. 
All  three  of  these  effects  can  be  taken  into  account  to  first  order  by  multiplying  by 
appropriate  scaling  factors. 

2.7.  Stability 

To  maintain  equilibrium,  the  electrical  and  mechanical  forces  must  be  equal.  By 
Equation  (3),  the  electrical  force  increases  quadratically  with  the  beam  deflection,  while 
the  mechanical  restoring  force  increases  linearly.  The  voltage  corresponding  to  the 
largest  stable  deflection,  known  as  the  pull-in  voltage,  can  be  found  by  equating  the 
mechanical  (Equation  (1))  and  electrical  (Equation  (3))  forces,  solving  for  the  voltage, 
and  finding  the  minimum  voltage  with  respect  to  the  distance  between  the  plate  and  the 
bottom  electrode,  resulting  in  the  following: 


(11) 

where  kejj is  the  total  effective  spring  constant  for  all  flexure  beams. 

The  steady-state  stability  of  the  canonical  FBMD  is  most  often  characterized  in 
terms  of  the  pull-in  point.  Going  a  step  further,  however,  it  can  be  shown  (by  equating 


23 


Equations  (1)  and  (3)  and  solving  for  g)  that,  for  simple  models  at  a  particular  voltage, 
exactly  two  steady-state  equilibrium  solutions  exist:  one  globally  stable  position  in  the 
gap  and  one  unstable  position.  The  latter  corresponds  to  a  deflection  past  one-third  of  the 
total  gap.  At  the  pull-in  point,  one  solution  exists,  beyond  which  all  positions  become 
unstable.  Physically,  this  phenomenon  is  a  result  of  the  inability  of  the  linear  mechanical 
spring  to  maintain  equilibrium  with  the  quadratic  electrostatic  force  throughout  the  gap; 
as  the  input  voltage  or  charge  increases,  causing  increased  beam  deflection  and  decreased 
electrode  separation,  both  the  electrical  attractive  force  and  mechanical  restoring  forces 
increase.  From  zero  deflection  to  the  pull-in  point,  the  mechanical  force  increases 
enough  to  maintain  steady-state  equilibrium  with  the  electrical  force,  but  after  this  point, 
the  mechanical  force  can  no  longer  increase  by  a  large  enough  amount  to  balance  the 
electrical  force  with  its  greater  rate  of  increase.  Having  lost  equilibrium,  the  top 
electrode  snaps  into  the  bottom  electrode  with  an  acceleration  proportional  to  the  net 
force,  resulting  in  the  eponymous  phenomenon. 

By  including  the  dynamics  of  the  system,  analytical  characterization  of  the  pull-in 
phenomenon  can  be  obtained.  Inertial  and  damping  forces  are  added  to  the  steady-state 
equation  as  in  the  following: 


-  (12) 

where  A,  B,  m,  k,  ,  and  V  are  the  mirror  plate  area,  damping  coefficient,  plate  mass, 
effective  spring  constant,  free-space  permittivity,  and  applied  voltage,  respectively.  The 
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variables  of  Equation  (12)  are  then  made  nondimensional  for  ease  of  analysis  through  the 
following  translations: 


(13) 


where  t  is  time,  to  produce 


—  (14) 

A  bifurcation  diagram  of  the  static,  non-dimensionalized  system  is  shown  in  Figure  6. 
The  dotted  line  corresponds  to  the  fold,  i.e.,  at  two-thirds  of  the  quiescent  gap  space,  go, 
only  one  C,  (non-dimensionalized  input  voltage)  exists.  The  shaded  gray  area  corresponds 
to  the  unphysical  region  where  the  deflection  is  greater  than  go.  As  asserted  above,  for 
each  voltage  input  C,  less  than  unity,  two  %  values  exist:  one  greater  and  one  less  than  one- 
third.  It  can  be  shown  that  the  former  is  an  unstable  equilibrium  and  the  latter  stable  [63]. 
The  existence  of  equilibria,  however  unstable,  beyond  a  third  of  the  gap  motivates  the  use 
of  closed-loop  control. 
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Figure  6:  Bifurcation  diagram  of  static,  non-dimensional  system 


2.8.  Kalman  Filter 

A  Kalman  Filter  (KF)  incorporates  sets  of  equations  describing  system  dynamics 
into  a  state-space  model,  which  it  combines  with  known  inputs,  initial  conditions,  and 
model  uncertainties  in  order  to  estimate  how  the  states — observable  or  not — and  their 
associated  covariances  change  over  time.  In  the  basic  form,  this  model  is  a  matrix  of 
linear,  first-order  differential  equations,  discretized  for  easy  use  with  digital  computers 
and  sample-and-hold  sensors.  Model  uncertainty  is  assumed  to  be  described  by  a 
Gaussian  distribution  with  zero  mean;  in  this  way,  only  the  first  two  statistical  moments 
need  monitoring.  Deviations  from  this  assumption,  such  as  biases  or  non-Gaussian 
processes,  can  be  modeled  by  augmenting  the  state  dynamics  model  with  noise  transfer 
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functions  driven  by  white  Gaussian  noise.  Estimates  of  the  states’  evolution,  propagated 
forward  in  time  using  the  dynamics  model,  are  refined  each  sample  period  by  reading  a 
measurement,  and  calculating  the  divergence  between  the  actual  and  expected  values  for 
that  measurement.  This  divergence,  known  as  the  residual  or  innovation  and  shown  in 
Figure  7,  provides  real-time  performance  data  of  the  KF;  modeling  inadequacies  will 
manifest  most  frequently  as  spikes,  biases,  or  divergence.  The  filter  gain  is  then  updated 
based  on  knowledge  of  uncertainties  in  the  measurement  process  and  dynamics  model  by 
taking  the  ratio  of  the  dynamics  covariance  and  the  residual  covariance.  The  product  of 
the  new  gain  and  residual  is  added  to  the  estimate,  generating  a  state  estimate  for  the  next 
sample  period,  which  serves  as  the  output  of  the  filter.  A  high  filter  gain  has  the  effect  of 
weighting  the  measurement  more  than  the  expected  value  in  the  KF  output,  while  a  low 
gain  has  the  opposite  effect. 

The  importance  of  the  residual  data  must  be  emphasized,  as  it  is  constitutes  the 
primary  means  by  which  the  adequacy  of  the  dynamics  model  may  be  quantitatively 
measured.  Real-world  system  deployments  do  not,  in  general,  enjoy  access  to  truth  data 
at  all,  much  less  in  real-time;  as  such,  absolute  estimation  errors  with  respect  to  truth  may 
not  be  quantified.  In  lieu  of  truth  data,  KF  designers  analyze  residual  data  for  indications 
of  model  inadequacies.  Residual  data  from  an  adequate  dynamics  model  will  feature  a 
zero  mean,  no  spikes,  and  root-mean-square  values  for  each  measured  state  equal  to  or 
less  than  system  tolerances  for  deviations.  Plots  of  residual  data  will  be  presented  in 
Chapter  5  for  the  current  system  and  analyzed  using  these  metrics. 
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Outputs  Should  Be 


Figure  7:  Block  diagram  illustrating  closed  loop  feedback  of 


the  KF 


2.9.  Extended  Kalman  Filter 

The  Extended  Kalman  Filter  (EKF)  expands  upon  the  operation  of  the  linear  KF 
by  taking  into  account  nonlinearities  inherent  in  the  dynamics  of  the  system,  without 
working  with  the  nonlinear  equations  directly.  Using  the  current  state  estimate  as  the 
operating  point,  the  EKF  employs  a  first-order  Taylor  polynomial  to  describe  the 
instantaneous  system  dynamics,  which  is  then  used,  along  with  covariance  and 
measurement  data  as  in  the  KF,  to  produce  a  new  state  output.  A  consequence  of  this  is 
that  the  matrices  required  to  generate  the  state  estimate  are  not  precomputable  as  they  are 
for  the  linear  KF,  thus  increasing  computation  time  considerably. 

The  advantage  of  this  scheme  is  greater  accuracy  in  estimating  strongly  nonlinear 
effects  than  a  KF  with  a  system  dynamics  model  linearized  about  a  static  operating  point; 
significant  errors  can  be  the  result  of  neglecting  strong  nonlinear  effects,  or  poor  choice 
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in  operating  point.  While  the  penalty  of  computation  is  undeniably  greater  than  that  of  a 
linear  KF,  it  remains  less  than  the  direct  use  of  nonlinear  equations,  with  acceptable 
accuracy  for  a  variety  of  applications. 

2.10.  Linear-Quadratic-Gaussian  (LQG)  Control 

The  LQG  controller  minimizes  a  quadratic  cost  function  assigned  to  some  or  all 
states  of  a  linear  or  linearized,  stochastic  system  that  is  disturbed  by  additive,  white 
Gaussian  noise.  Continuous  or  discrete  measurements,  functionally  related  to  some  or  all 
states,  are  assumed  to  be  corrupted  by  additive,  white  Gaussian  noise  of  zero  mean  and 
fed  back  to  the  controller.  Although  assumed  to  be  not  applicable  here,  non-zero-mean 
noise  sources  may  be  modeled  by  an  appropriate  shaping  filter  driven  by  zero-mean, 
white  Gaussian  noise.  For  the  purposes  of  this  work,  discrete-time  sensor  data  is 
assumed  to  be  available  as  the  analog-to-digital  conversion  of  sampled,  continuous-time 
data. 

The  cost  function  is  defined  as  the  first  moment  of  the  sum  of  the  states  and 
control  inputs  over  all  time  and  the  desired  final  state,  each  squared  and  weighted 
according  to  the  constraints  of  the  design  and  the  tractability  of  the  problem. 
Equivalently,  states  can  be  replaced  by  expressions  quantifying  deviation  from  a 
reference  trajectory.  In  either  case,  control  inputs  are  found  such  as  to  minimize  the  cost 
function,  with  a  conceptual  “control  energy”  commensurate  with  their  respective  weights. 
That  is  to  say,  the  larger  the  weight,  the  larger  the  control  effort  is  expended  in  driving 
the  state  or  deviation  to  zero,  or  the  smaller  the  control  input.  The  former  is  useful  in 
accurate  trajectory  tracking  applications,  for  example,  while  the  latter  might  conserve 
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finite  available  energy  or  smooth  the  controller  response. 

Minimization  of  the  cost  function  is  accomplished  by  solving  a  trio  of  Riccati 
differential  equations,  two  running  forward  in  time  and  one  running  backward.  Those 
running  forward  model  the  state  and  covariance  dynamics  and  constitute  the  estimator 
problem,  as  solved  by  the  KF  and  EKF  described  above  and  shown  in  Figure  7  as 
“Optimal  State  Estimate.”  The  backward-running  equation  solves  for  the  feedback  gain 
matrix,  the  product  of  which,  with  the  state  estimate  (i.e.,  the  output  of  the  KF), 
determines  the  gains  in  the  general  control  law  to  be  applied  to  the  system,  as  well  as  the 
“Model  of  System  Dynamics”  box  in  Figure  7.  This  input  is  submitted  to  a  zero-order 
hold  and  applied  at  the  start  of  the  next  sample  period.  The  control  law  for  this  work  is 
derived  and  specified  in  Chapter  4. 
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3.  Modeling 


3.1.  Chapter  Overview 

An  accurate  model  is  essential  to  the  function  of  an  observer  such  as  the  Kalman 
Filter.  The  observer  estimates  system  states — such  as  velocity  or  accumulated  charge — 
and  disturbances  using  knowledge  from  the  model  and  about  the  passage  of  time. 
Typically,  as  it  is  here,  the  model  is  a  state-space  representation  of  the  relevant  equations 
of  motion. 

This  chapter  is  devoted  to  the  development  of  such  a  model  for  the  micromirror 
described  in  Chapter  2  by  means  of  simulation  and  analysis.  Various  nonlinearities 
particular  to  the  FBMD  under  consideration  were  introduced  in  Chapter  2  and  will  be 
quantified  here  through  simulation.  First,  a  pull-in  voltage  study  is  perfonned  in 
CoventorWare,  from  which  parasitic  capacitance  and  effective  spring  constant  are 
derived.  Next,  squeeze-film  damping  is  modeled  in  COMSOL  Multiphysics.  These 
effects  are  combined  into  a  governing  equation  of  motion,  which  is  solved  in  MATLAB 
in  order  to  characterize  nonlinear  transients. 

3.2.  Pull-In  Voltage 

CoventorWare  contains  an  algorithm  to  detect  pull-in  conditions  during  its 
simulation  of  electromechanically  defonned  geometries.  When  engaged,  the  algorithm 
checks  the  calculated  displacement  for  divergence,  indicating  a  pull-in  condition. 
CoventorWare  then  simulates  an  input  equal  to  the  mean  of  the  last  converging  and 
diverging  inputs,  iterating  until  either  the  simulation  takes  a  maximum  number  of  steps, 
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or  the  result  is  less  than  a  set  tolerance — that  is,  a  change  in  calculated  deflection.  The 
result  of  this  iteration  constitutes  the  last  point  of  stable  equilibrium  before  pull-in.  The 
pull-in  voltage  for  the  present  study’s  nominal  FBMD  is  calculated  in  this  manner. 

A  three-dimensional  model  of  the  FBMD  is  created  first  by  drawing  its  layout  in 
Tanner  L-Edit.  This  layout  file  specifies,  based  on  user  input,  the  length  and  width  of  all 
desired  layers  in  the  Poly  MUMPS  process.  The  file  may  be  imported  into  CoventorWare 
and,  when  combined  with  the  appropriate  PolyMUMPS  fabrication  process  definition, 
which  specifies  the  process  described  in  Chapter  2  in  a  standardized  way,  CoventorWare 
creates  a  three-dimensional  model,  complete  with  material  definitions  and  layer 
thicknesses  (Figure  8).  Note  the  conformal  nature  of  the  process  reflected  in  the  shape  of 
each  of  the  flexures  via  the  rectangular  “bump”  near  the  supports. 


Figure  8:  Three-dimensional  meshed  CoventorWare  model  of  the  FBMD  with  an 
input  voltage  of  8  V. 
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The  pull-in  study  begins  by  specifying  mechanical  and  electrical  boundary 
conditions:  the  polyO  layer  is  immobile  and  connected  to  the  electrical  ground;  the  polyl 
flexure  supports  are  immobile;  no  external  mechanical  load  is  applied  to  the  mirror  plate 
or  flexures;  no  mechanical  contact  between  surfaces;  and  polyl  is  connected  to  the  input 
voltage.  As  such,  the  polyl  flexures  and  mirror  plate  are  free  to  respond  to  electrical  and 
mechanical  forces  and  move  in  space.  With  boundary  conditions  in  place,  an  input 
voltage  trajectory  is  set  in  CoSolveEM  to  be  applied  to  the  FBMD’s  input — in  this  case,  a 
voltage  ramp  from  0  V  to  25  V  in  0.5  V  steps.  This  voltage  range  was  detennined  by 
solving  Equation  (11)  with  the  FBMD  dimensions  specified  in  Chapter  2  and  the 
parameters  measured  by  MUMPS  in  Run  77  (composite  Young’s  Modulus  E  of  131 
GPa),  yielding  an  analytically  expected  pull-in  voltage  of  18.35  V.  Equation  (11)  does 
not  take  into  account  any  of  the  non-idealities  described  in  Chapter  2,  so  the  simulated 
voltage  range  was  well  in  excess  of  the  analytical  pull-in  point  to  ensure  pull-in  was 
captured.  Also  note  that  the  CoventorWare  simulation  assumes  steady-state  conditions. 

The  simulation  is  executed,  and  the  pull-in  voltage  is  found  to  be  21.44  V  (Figure 
9).  For  each  input  voltage  in  the  trajectory  up  to  pull-in,  the  corresponding  mechanical 
force,  electrical  force,  total  charge  accumulated  in  polyl,  total  capacitance,  and  deflection 
are  calculated.  These  results  show  strongly  nonlinear  mechanical  force  versus  voltage,  as 
well  as  nonlinear  charge  versus  voltage,  indicating  a  nonlinear  parasitic  capacitance. 
Each  of  these  will  be  discussed  in  the  following  sections. 
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Figure  9:  Deflection  and  Total  Charge  versus  input  voltage;  pull-in  occurs  at 
21.44  V 

3.3.  Effective  Spring  Constant 

In  Chapter  2,  a  simple  analysis  based  on  Euler-Bemoulli  equations  for  beam 
bending,  shown  in  Equation  (10),  was  presented  to  be  used  as  a  lumped  model  parameter 
in  a  mass-spring-damper  construct.  This  equation  yields  a  value  of  12.576  N/m  for  an 
effective  linear  spring  constant  using  the  same  Young’s  Modulus  as  above  of  131  GPa. 
This  analytical  value  is  now  compared  against  a  full  three-dimensional  simulation  in 
CoventorWare,  which  includes  the  rigid  support  stacks,  hinges,  and  residual  and  internal 
stresses.  To  enable  a  direct  comparison,  the  simulated  mechanical  force  is  plotted  against 
geometrical  deflection  from  rest  height,  on  which  a  linear  regression  is  performed, 
forcing  the  zeroth  order  term  to  zero  (Figure  10).  The  slope  of  this  line,  corresponding  to 
the  simulated  effective  spring  constant,  is  found  to  be  15.056  N/m,  with  a  correlation 
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Figure  10:  Simulated  mechanical  restoring  force  versus  deflection; 
effective  k=15.059  N/m 

coefficient  of  0.9999.  The  latter  quantity  belies  indication  of  a  high  degree  of  linearity  of 
the  line,  as  the  regression  calculation  is  dominated  by  a  preponderance  of  points  at  small 
deflection,  tending  to  be  more  linear.  Divergence  from  the  line  increases  with  deflection, 
but  in  the  domain  of  interest,  the  effect  remains  small  enough  to  be  neglected. 

Several  non-ideal  mechanisms  contribute  to  this  discrepancy,  one  of  which  is  the 
restoring  force  of  the  hinge  connecting  the  flexure  beam  to  the  mirror  plate.  This  hinge  is 
8  |am  in  the  direction  parallel  to  the  length  of  the  flexure  (w),  and  3  qm  perpendicular  (/), 
creating  a  3  qm  space  between  flexure  and  mirror.  As  the  flexure  deflects  from  the 
nominal  height  and  the  mirror  descends,  a  torque  is  introduced  in  the  hinge,  which  adds 
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to  the  mechanical  potential  energy  of  the  flexure  beams.  This  deflection  is  never  more 
than  2%  of  the  length  L  of  the  flexure,  thus  in  the  limit  of  small  angles,  the  ratio  of  the 
vertical  displacement  to  the  flexure  length  approximates  the  angle  over  which  the  torque 
is  applied.  Assuming  the  force  inducing  the  torque  is  applied  perpendicularly  to  the  lever 
ann  (in  this  case,  the  flexure  beam),  the  equation  from  [19]  is  modified: 


(15) 

where  w  is  the  width  of  the  hinge  as  above;  t  the  thickness  (2  pm);  E  is  Young’s 
Modulus;  /  the  length  of  the  hinge  as  above;  d  the  maximum  distance  of  the  flexure 
deflection;  and  //  is  Poisson’s  Ratio  for  polysilicon  (0.22  [67]).  This  equation  models  the 
hinge  as  another  linear  spring  in  series  with  the  flexures,  thereby  adding  a  small,  but  non¬ 
zero,  value  of  0.013  N/m  to  the  effective  spring  constant — the  lumped  model 
parameter — of  the  system. 

Another  stiffening  effect  is  that  of  two-dimensional  strain,  commonly  referred  to 
as  the  plate  effect.  As  the  ratio  of  flexure  width  to  length  increases,  so  too  does  the 
curling  of  the  beam  in  the  axis  of  width,  around  the  axis  of  length,  creating  what 
resembles  a  valley  or  half-pipe,  with  a  curvature  opposite  in  sign  to  that  of  the  flexure 
deflection.  The  net  result  of  this  transverse  curling  is  in  effect  to  stiffen  the  flexure,  and  is 
usually  modeled  by  use  of  a  biaxial  modulus  in  place  of  Young’s  modulus  in  the 
analytical  beam-bending  equations.  This  biaxial  modulus  is  calculated  by  dividing  the 
Young’s  modulus  by  difference  of  unity  and  Poisson’s  ratio  of  the  material,  and,  since  the 
latter  is  always  less  than  one,  the  biaxial  modulus  is  always  greater  than  the  Young’s 
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modulus.  The  Euler-Bernoulli  equations  can  be  re-derived  to  remove  this  dependence  on 
Poisson’s  ratio,  allowing  direct  insight  into  the  effect  of  the  width-to-length  ratio.  As  the 
width-to-length  ratio  for  the  system  currently  studied  is  0.03,  it  is  found  [70]  that  the 
plate  effect  is  likely  not  a  source  of  deviation. 

Another  possible  source  of  effective  spring  constant  deviation  is  the  increase  of 
structural  compliance  if  the  design  features  built-up,  pillar-like  supports,  such  as  the 
flexure  beam  anchors  described  for  the  current  FMBD  in  Chapter  2.  These  supports  are 
able  to  rotate  in  the  presence  of  external  moments  induced  by  residual  stress.  It  is  found 
that  this  can  be  analytically  modeled  [70]  as  a  small  (relative  to  that  designed)  increase  or 
decrease  in  mechanical  flexure  length.  That  is  to  say,  a  beam  attached  to  an  anchor  that 
rotates  behaves  similarly  to  a  slightly  longer  beam.  Finite  element  modeling  results  and 
analytical  estimations  in  [70]  for  maximum  cantilever  deflection  versus  several  drawn 
mask  lengths  were  plotted.  It  was  found  that  the  addition  of  5.85  micron  to  the  cantilever 
design  length  in  the  analytical  equations  fits  the  FEM  data  well.  For  this  work, 
subtracting  5.85  micron  to  the  nominal  design  length  of  100  micron  yields,  via  Equation 
(10),  an  expected  mechanical  spring  constant  of  15.069,  within  10%  of  that  simulated  in 
CoventorWare.  This  structural  compliance,  combined  with  the  addition  of  torsion  spring 
action  as  described  above,  are  then  the  dominant  effects  in  deviating  from  the  first-order 
estimation  in  Equation  (10). 

3.4.  Parasitic  Capacitance 

This  section  serves  to  summarize  modeling  of  the  capacitance  of  micromirrors  as 
simulated  in  CoventorWare.  To  start,  the  micromirror  was  simulated  for  input  voltages 
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ranging  from  zero  to  pull-in,  the  latter  being  estimated  by  Coventor’s  pull-in  detection 
algorithm.  For  each  voltage,  a  total  accumulated  charge,  deflection  distance,  and 
capacitance  was  recorded  and  exported.  Excel  was  then  used  to  tabulate  the  results  for 
comparison  to  theory  and  fitting. 

As  an  initial  sanity  check,  the  exported  capacitance  values  were  plotted  against 
the  ratio  of  accumulated  charge  and  voltage  for  each  input  voltage,  and  found  to  match 
exactly.  The  former  was  then  plotted  against  the  familiar  equation  for  the  capacitance  of 
infinitely  long  parallel  conductors,  employing  the  mirror  plate  area  and  exported 
deflection  distances.  Stark  differences  between  the  exported  capacitances  and  theory 
were  seen,  both  magnitude  and  rate  of  increase  near  the  pull-in  voltage. 

A  lumped  parameter  model  was  sought  to  include  these  nonlinearities,  likely  the 
results  from  neglecting  the  fringe  field  effects  and  capacitive  contributions  from  the  four 
flexure  beams.  This  model  was  to  take  the  fonn  of  the  original,  first-order  equation  for 
parallel  plates,  but  using  an  effective  area  that  would  account  for  the  aforementioned 
nonlinearities.  Both  this  effective  area  model  and  that  of  exported  capacitance  were 
strongly  nonlinear  near  the  pull-in  voltage  and  converge  with  increasing  input  voltage; 
the  maximum  divergence  of  3.945  fF  (roughly  one  order  of  magnitude  smaller)  occurs  at 
3.50  V.  A  misguided  attempt  to  fit  an  analytical  equation  to  this  divergence  curve 
resulted  in  a  power  law  with  a  correlation  coefficient  of  0.997.  This  approach  was 
abandoned  due  to  complexity.  Instead,  a  linear  regression  against  the  simulated 
capacitance  versus  gap  data  was  perfonned  and  found  to  model  more  closely  the 
difference  between  theory  and  observed  capacitance  over  the  range  of  gap  values  (Figure 
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Figure  11:  Simulated  post-process  capacitance;  ratio  of  simulated  total  charge  and 
voltage;  differential  charge  over  voltage  ratio;  Leus's  fringe  field  corrected  capacitance 
[71];  and  basic  parallel  plate  capacitance  versus  input  voltage 

11).  The  correlation  coefficient  for  this  linear  fit  was  0.9999.  Separately,  Microsoft 
Excel’s  linear  least-squares  regression  function,  LINEST,  was  used  to  fit  a  linear 
equation  directly  (in  contrast  to  the  difference  curve  above)  to  the  exported  capacitance 
versus  the  inverse  of  the  exported  gap,  with  y-intercept  forced  to  zero.  Dividing  the  fitted 
slope  by  the  free-space  permittivity  yielded  the  effective  area:  1.07x10'  m  . 

The  contributions  of  fringe  fields  and  the  flexures  were  investigated  in  an  attempt 
to  explain  observed  differences  between  theory  and  simulation.  First,  the  flexures  deflect 
according  to  the  Euler-Bemoulli  equation  (Equation  (9))  presented  in  Chapter  2,  and  can 


39 


be  considered  as  capacitors  with  one  planar  plate  and  one  curved  plate.  Since  the  inverse 
of  this  equation  is  difficult  to  integrate  analytically,  an  approach  emulating  a  Riemann 
sum  is  taken  to  model  the  flexures  as  a  series  of  parallel-plate  capacitors,  connected  in 
parallel,  with  gap  equal  to  the  amount  of  deflection  at  that  distance  away  from  the  fixed 
support.  As  this  is  a  first-order  attempt  to  quantify  the  order  of  magnitude  of  flexure 
capacitance,  fringe  fields,  non-uniformities  introduced  by  fabrication,  and  the  fact  that  the 
ground  plane  is  not  directly  below  the  flexure  are  ignored.  The  analysis  proceeds  by 
dividing  the  flexure  into  ten  equal  parts,  calculating  the  deflection  at  the  midpoint  of  each 
part  to  be  used  as  the  gap  between  plates,  inserting  this  gap  into  the  capacitance  equation 
for  parallel  plates,  and  adding  the  contributions  for  the  ten  sections  and  four  flexures. 
The  resulting  contribution  is  two  orders  of  magnitude  smaller  than  that  of  the  mirror 
plate,  and  therefore  may  be  neglected. 

Fringe  fields  are  modeled  by  previous  work  [71]  in  an  analytical  equation  that 
accounts  for  non-zero  plate  thickness,  side,  and  edge  effects.  The  equation  itself  is  a 
modified  version  of  earlier  work,  the  result  of  which  is  a  Schwartz-Christoffel  conformal 
mapping  transformation  describing  the  electric  field  with  parallel-plate  boundary 
conditions.  What  results  underestimates  the  simulated  capacitance  at  low  voltages  (larger 
gaps)  and  overestimates  at  high  voltages  (smaller  gaps),  but  only  by  a  maximum  of  three 
percent,  outperfonning  both  the  basic  parallel  plate  theory  and  the  parallel  plate  equation 
linearly  fitted  effective  area.  This  performance  supports  the  need  to  consider  fringe  fields 
in  an  accurate  treatment  of  the  total  capacitance. 
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The  differential  definition  of  capacitance,  i.e.,  the  ratio  of  the  change  in  charge 
and  the  change  in  input  voltage,  is  necessary  as  the  basic  theory  for  parallel  plates 
assumes  a  static  and  uniform  gap  between  the  plates;  since  this  gap  has  a  voltage 
dependence  (indeed,  the  modus  operandi  of  electrostatic  actuation),  the  differential 
definition  of  capacitance  is  more  appropriate.  To  generate  this  differential  data  for 
comparison,  a  linear  voltage  trajectory  was  again  applied  to  the  device,  offset  from  the 
original  trajectory  by  0.1  V,  which  enables  a  quasi-static  differential  analysis.  As  a 
second  sanity  check,  this  differential  model  and  the  original  capacitance  data  were  both 
plotted  against  input  voltage  and  found  to  be  strikingly  dissimilar,  providing  some 
evidence  for  the  inadequacy  of  CoventorWare’s  use  of  the  steady-state  definition  of 
capacitance,  which  is  the  total  accumulated  charge  divided  by  potential.  CoventorWare 
obtains  each  of  the  values  in  this  ratio  by  first  inferring  a  surface  charge  on  conductive 
surfaces  based  on  changes  in  the  electric  displacement  vector,  and  then  integrating  over 
the  total  surface  area  of  the  structure,  yielding  a  total  charge.  The  latter  is  then  divided 
by  input  voltage  and  exported  as  a  total  capacitance.  This  can  be  seen  by  simply  dividing 
the  exported  total  charge  by  the  known  input  voltage,  which  exactly  follows  the  exported 
capacitance  curve.  The  trouble  with  this  approach  is  the  relationship  between  the  charge 
and  the  potential.  The  latter  is  defined,  to  first-order,  as  the  strength  of  the  electric  field, 
E,  set  up  between  charges  on  the  mirror  plate  and  ground  plane,  multiplied  by  the 
distance  between  the  charges,  d,  i.e.,  V=Ed.  Since  this  device  is  electrostatically 
actuated,  the  accumulated  charge  interacts  with  this  distance  by  design;  in  particular,  a 
feedback  loop  is  set  up  between  the  distance  between  the  plates  and  the  amount  and 
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location  of  charge  on  each  plate.  To  assume  a  constant  distance  is  a  mistake,  the 
consequence  of  which  is  to  ignore  fundamental  dynamics  of  the  system.  As  such,  the 
capacitance  values  given  by  CoventorWare  will  be  considered  suspect. 


3.5.  Fluid  Damping 

In  order  to  characterize  any  film  damping  the  micromirror  might  experience,  a 
three-dimensional,  one-quarter  model  of  the  mirror  plate — with  the  remainder  modeled 
by  numerical  continuity — is  built  in  COMSOL  Multiphysics  simulator.  A  perfectly 
diffuse  tangential  momentum  accommodation  coefficient  for  the  fluid  is  also  assumed; 
this  is  reasonable  based  on  work  [69]  on  gas  flows  over  pure  silicon,  which  demonstrated 
diffuse  accommodation  despite  orders  of  magnitude  smoother  surfaces.  Slip  flow  was 
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Figure  12:  Damping  force  (N)  and  mirror  deflection  (m)  vs.  time  (s) 
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accounted  for  in  the  simulation  as  well,  since  the  system  has  a  Knudsen  number  of  0.034. 
A  constant  potential  is  specified  across  the  mirror  and  ground  plates,  and  the  simulation 
is  conducted  to  record  the  transient  response  of  all  applicable  physical  forces: 
mechanical,  electrical,  and  fluid  (Figure  12).  The  velocity  of  the  plate  at  each  sample 
period  is  calculated  numerically  in  post-processing  using  a  simple  difference  function  in 
MATLAB.  The  COMSOL-calculated  damping  force  is  then  divided  at  each  sample 
period  by  the  plate  velocity  to  generate  an  estimate  of  the  damping  coefficient,  B ,  at 
1.1 885x1  O'5  Ns/m.  As  the  plate  settles  to  equilibrium,  the  velocity  oscillates  around  zero, 
causing  the  estimate  of  B  to  vary  severely;  as  such,  only  a  small  subset  of  samples  around 
the  maximum  velocity  are  used  to  calculate  B ,  as  they  displayed  the  most  stable  velocity 
values. 

3.6.  Transient  Analysis 

The  above  analyses  quantify  various  physical  characteristics  required  to  describe 
the  FBMD  with  a  system-level  model,  which  will  be  directly  employed  by  the  Kalman 
Filter  (developed  in  Chapter  4)  to  estimate  unobservable  states — specifically,  and  for 
reasons  outlined  in  Chapter  1,  the  deflection  of  the  mirror  plate.  As  deflection  and 
potential  are  interdependent  upon  each  other  and  demonstrate  strong,  nonlinear  transients, 
these  characteristics  are  simulated  together  to  verify  the  claims  of  stability  in  Section  2.7 
that  were  based  on  simple,  static  models.  The  results  are  shown  in  Figure  13.  Despite 
beam  and  fluid  damping  transients,  delays  for  charge  accumulation  (albeit  minute),  and 
all  the  nonlinearities  described  above,  the  results  indicate  stable  equilibria  are  possible  up 
to  just  over  18.3  Volts.  This  value  is  decidedly  close  to  the  analytically  derived  value  of 
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18.4  Volts  in  Chapter  2,  but  assumes  perfect  polysilicon  edges  and  does  not  account  for 
the  “bump”  in  each  of  the  flexure  beams  as  CoventorWare  does.  Beyond  the  sanity 
check,  the  transients  also  indicate  that  stability  is  reached  within  15  to  20  microseconds, 
and,  depending  on  the  magnitude  of  the  input,  instability  begins  running  away  before  40 
microseconds.  This  implies  the  sample  rate  of  the  controller  must  be  50  kHz  to  meet  the 
Nyquist  criterion,  and  250  kHz  for  a  margin  often  samples  on  the  transient. 
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Figure  13:  Mirror  plate  deflection  (m)  vs.  time  (s)  via  COMSOL.  Input  voltage  swept  from 
zero  to  18.34  V  (top,  diverging  trajectory). 
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4.  State  Estimator  and  Controller  Design 


4.1.  Chapter  Overview 

This  chapter  describes  the  designs  for  the  observer  and  controller  algorithms,  to 
include  assumptions  made  in  the  measurement  process,  and  concludes  with  an  overall 
algorithm  description.  It  is  split  into  three  interdependent  parts:  the  first  for  the  observer, 
the  second  for  the  controller,  and  the  third  for  the  algorithm  description.  The  first  part 
concludes  with  process  equations  used  for  the  observer.  The  second  section  concludes 
with  the  control  law  equations.  The  last  part  steps  through  the  MATLAB  scripts  that 
execute  the  first  two  parts. 

4.2.  State  Estimator  Design 

The  controller  requires  an  accurate  estimate  of  the  controlled  variable  to  be 
effective.  The  observer,  as  such,  is  the  lynchpin  to  effective  control,  as  it  provides  an 
estimate  of  the  controlled  variable  in  spite  of  indirect  observation  and  noise.  This  section 
describes  the  design  of  a  Kalman  Filter  to  be  applied  to  the  FBMD  as  analyzed  in 
previous  chapters. 

4.2.1.  Measurement  Design  and  Assumptions 

Using  parallel-plate  capacitor  charge  accumulation  and  mirror  plate  dynamics 
relationships  developed  in  Chapters  2  and  3,  the  observer  infers  a  mirror  plate  deflection 
from  measured  changes  in  voltage  across  the  mirror.  The  proposed  measurement  system 
is  modeled  as  an  ideal  voltmeter,  to  include  an  internal  resistance  of  1  Mf2.  This  design 
is  motivated  by  the  ease  of  implementation  in  a  lab  for  empirical  verification,  and  is 
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demonstrative  of  a  general  system  of  electrical  measurement.  To  obtain  the  voltage 
across  the  FBMD,  the  mirror  is  driven  by  the  series  connection  of  an  ideal  voltmeter  (in 
particular  an  independent  voltage  source  and  an  internal  resistor)  and  the  FBMD.  The 
Norton  equivalent  of  this  circuit  is  employed  (Figure  14),  making  the  calculation  of  the 
amount  of  current  going  to  the  FBMD  a  simple  application  of  Kirchoff  s  Current  Law 
(KCL):  the  current  to  the  FBMD  is  the  difference  of  the  current  source  (control  signal) 
and  the  current  through  the  internal  resistor.  Since  the  current  through  the  FBMD  is  the 
time-differential  of  the  voltage  across  it  scaled  by  the  effective  capacitance  at  that  instant 
(Equation  (16)),  the  observable  state  variable  is  obtained  as  the  integral  of  the  difference 
of  the  control  signal  and  the  resistor  current,  divided  by  the  effective  capacitance  at  that 
moment: 


— -  (16) 

-  -  (17) 

Modeling  the  FBMD  as  a  voltage-controlled  capacitance  in  this  way  assumes  that  the 


Figure  14:  Norton-equivalent  driving 
circuit.  FBMD  shown  here  as  a  capacitor. 
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capacitance  is  linear,  allowing  the  time  derivative  of  the  capacitance  to  go  to  zero. 

As  this  observer  will  be  implemented  digitally,  the  voltmeter  will  sample  the 
voltage  across  the  FBMD  and  hold  for  one  sample  period.  This  process  is  repeated  every 
50  nanoseconds,  modeling  a  voltmeter  capable  of  20  MHz  operation.  The  voltmeter  is 
subject  to  thermal  drift  and  imprecisely  known  internal  resistance,  but  performs  no 
filtering  or  buffering  of  more  than  one  sample.  Other  sources  of  noise  (shot  noise,  cross¬ 
talk,  etc.)  may  be  added  with  appropriate  models;  for  generality,  the  noise  in  the  filter 
algorithm  is  modeled  as  additive,  white,  and  Gaussian,  with  a  magnitude  of  0.01  and 
mean  of  zero. 

In  this  work,  measurements  are  generated  in  software  by  the  MATLAB  script 
get_meas().  Inputs  to  the  script  include  absolute  start  and  stop  times  (the  products  of  the 
sample  period  and  the  previous  and  current  sample  numbers),  the  previous 
measurements,  and  the  previous  control  signal,  the  script  outputs  the  current  sample 
period’s  measurement  of  voltage  across  the  FBMD,  as  well  as  the  gap  measurement  for 
troubleshooting.  The  script  calculates  the  measurements  using  MATLAB ’s  ode45  script 
and  the  nonlinear  functions  developed  in  Chapter  3,  with  the  boundary  conditions 
specified  by  the  time  period  taken  as  input  by  get_meas(),  initial  state,  and  control  signal. 
This  implementation  allows  for  changing  control  signals  in  real-time,  at  the  severe  cost  of 
calling  ode45()  for  each  sample  period;  pre-calculation  of  the  control  history  would  speed 
up  the  script  significantly  by  enabling  the  measurement  history  to  be  pre-calculated  and 
called  as  a  look-up  table.  This,  however,  would  not  be  representative  of  a  real-time 
controller  implementation. 
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4.2.2.  Propagation  Design 

With  the  measurement  system  described  in  the  previous  section,  and  equations  of 
motion  analyzed  in  Chapter  3,  the  following  state  dynamics  equations  describe  using 
nonlinear,  first-order  partial  differential  equations  the  system  states  as  potential  across 
mirror  plates,  mirror  deflection,  and  mirror  velocity: 


(18) 

The  state  noise  covariance,  quantified  by  the  matrix  Q,  indicates  the  uncertainty  in  the 
dynamics  equations.  Real-world  disturbances,  such  as  drift,  wander,  and  non- white 
noise,  would  require  additional  states  modeled  by  shaping  filters  driven  by  white 
Gaussian  noise.  In  this  work,  no  assumptions  may  be  made  about  the  deployed 
environment;  rather,  pseudo-noise  is  added  the  dynamics  equations  in  an  attempt  to 
account  for  unmodeled  effects.  Q  can  be  adjusted  to  accommodate  unexpected 
disturbances  and  system  model  inadequacies,  such  as  the  consequences  of  linearizing 
dynamics  equations  of  higher  polynomial  orders.  The  loss  of  the  effects  of  higher  order 
tenns  may  be  conceptually  considered  by  the  algorithm  as  an  unmodeled  disturbance  and 
negotiated  in  the  same  way,  by  increasing  the  state  noise  covariance  Q.  This  increase  has 
the  effect  of  adding  weight  to  the  measured  state  value  over  the  propagation  calculation 
in  the  observer  algorithm,  desirable  during  the  initial  transient  or  large  disturbances.  Too 
much  weight  on  measurements  vice  calculated  estimates  can  be  detrimental  to  a 
detectable  system.  The  measurement  improves  the  estimate  of  the  controlled  variable 
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(FBMD  deflection)  using  a  voltage  in  conjunction  with  modeled  dynamics;  too  little 
weight  on  the  latter  creates  systemic  errors  that  the  Kalman  Filter  has  no  means  to 
remove.  By  trial  and  error,  it  was  found  that  values  of  lxl O'4,  lxl O'9,  and  zero  along  the 
diagonal  of  Q  (roughly  0.1%  of  each  expected  state  magnitude)  provided  a  good 
compromise. 

Pursuant  to  the  assumption  that  small  deviations  from  system  equilibrium  can  be 
treated  as  linear,  the  dynamics  equations  are  linearized  before  discretization.  Equations 
for  equilibrium  state  and  control  values  are  found  as  functions  of  the  desired  setpoint  by 
setting  the  dynamics  equations  to  zero  and  solving  for  each  state  and  for  the  control 
signal: 
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The  dynamics  equations  are  linearized  by  calculating  the  Jacobian  and  setting  the  states 
and  control  signal  to  their  respective  equilibrium  values  per  equations  previously 
described: 
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(21) 


(22) 

Q  is  as  described  in  the  previous  section,  and  G  is  set  to  the  identity  matrix.  These 
matrices,  F,  B,  G,  and  Q,  along  with  the  sample  period  At,  are  input  to  the 
equiv_discrete()  script,  which  employs  the  MATLAB  c2d()  function  to  generate  a 
discretized  state  dynamics  matrix,  O,  control  matrix  Bd,  and  state  dynamics  noise  matrix 
Qd,  from  respective  continuous  equations,  each  of  which  are  used  to  propagate  the 
estimate  of  state  in  the  Kalman  Filter  from  one  sample  period  to  the  next. 

Initial  conditions  for  each  filter  simulation  are  assumed  to  be  at  quiescence,  in  this 
case,  zero  for  all  states  and  control.  In  this  case,  MATLAB  finds  difficulty  numerically 
solving  the  differential  equations  starting  at  pure  zero,  thus  “near”  zero  values  of  1x10' 15 
are  used  instead.  The  state  covariance  matrix  P  has  an  initial  value  of  10  along  the 
diagonal.  This  value  was  found  by  trial  and  error  to  have  an  appropriately  small  transient 
behavior. 

4.2.3.  Truth  Generation  Algorithm 

Filter  perfonnance  may  be  judged  by  analysis  of  residuals  and  error  history. 
Residuals  are  calculated  by  the  difference  of  the  measured  values  and  post-propagation 
state  estimates  mapped  to  measurement  space.  As  described  above,  the  measured  values 
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are  not  empirical  data,  but  rather  are  simulated  by  periodic  solutions  to  nonlinear, 
ordinary  differential  equations  (ODE) — the  same  equations  from  which  the  filter’s 
dynamics  model  are  linearized.  Although  residual  analysis  is  typically  very  useful,  it 
merely  serves  in  this  case  to  indicate  how  well  the  linearized  filter  predicts  the  output  of  a 
more  robust  ODE  solver.  Error  history  analysis  is  essential  in  this  case,  and  requires  a 
dataset  to  act  as  truth  for  comparison. 

The  FEM  analysis  as  executed  in  COMSOL  serves  as  this  dataset.  As  described 
in  Chapter  3,  the  FBMD  system  analysis  consisted  of  setting  each  plate  of  the  FBMD  at  a 
constant  potential  and  recording  the  time  history  of  charge  accumulation  and  mirror 
deflection.  By  contrast,  the  filter  applies  a  known  (not  necessarily  constant)  current  and 
measures  the  change  in  voltage  over  time.  Furthermore,  the  FEM  and  filter  time  vectors 
do  not  overlap.  This  disparity  is  overcome  by  considering  the  FEM  dataset  as  a  look-up 
table,  with  geterr()  as  a  look-up  script,  rather  than  a  directly  analogous  deflection 
trajectory  to  be  compared  point-for-point  with  the  observer’s  results.  The  raw  FEM  data 
is  first  organized  into  structured  arrays  by  electrical  potential,  while  removing  redundant 
data  necessary  for  FEM  accuracy.  The  MATLAB  timeseries  command  is  employed  to 
create,  for  each  FEM- simulated  electrical  potential,  a  time  trajectory  of  maximum  mirror 
deflection,  with  which  the  resample  command  is  used  to  match  observer  and  FEM  time 
vectors.  Resample  introduces  some  small  amount  of  error  in  the  data,  as  it  estimates  data 
by  linearly  interpolating  between  times  in  the  FEM  data  to  align  with  those  in  the 
observer.  Once  the  FEM  data  are  resampled  along  the  observer’s  time  vector,  geterr() 
calculates  the  total  charge  accumulated  on  the  FBMD  plates  for  each  time  sample,  then 
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for  the  same  voltage  finds  the  charge  and  deflection  as  determined  by  the  FEM.  The 
difference  of  FEM  and  observer  charge  and  deflection  values  for  each  time  is  output  as 
the  error  array. 

4.3.  Controller  Design 

The  central  goal  of  this  work  is  to  mitigate  mirror  snap-in  and  demonstrate 
controlled  deflection  of  the  FBMD  throughout  the  gap.  To  this  end,  the  detenninistic 
control  law  sought  will  achieve  and  maintain  (requiring  type-1  for  modeling  errors)  a 
desired  non-zero  setpoint,  and  reject  disturbances.  Typically,  one  or  more  system-level 
performance  requirements  would  be  levied  to  guide  the  control  design.  For  example, 
consider  a  free-space  optical  communications  system  employing  an  array  of  FBMDs  for 
wavefront  error  correction;  characteristics  such  as  deflection  accuracy  of  +/-  1  nm,  1  ns 
response  time  to  steady  state,  and  10%  overshoot  may  be  necessary  for  an  acceptable  bit 
error  rate,  enabling  mission  success.  In  this  work,  design  choices  were  made  primarily  to 
mitigate  snap-in,  but  system-level  considerations  will  be  indicated  throughout  the  design 
description. 

4.3.1.  Pseudo  Rate  Au 

One  means  to  effect  integral  control  is  to  cost  not  only  differences  in  control 
signal  from  the  nominal  uo  calculated  above,  but  also  differences  between  samples  in 
time.  The  state  vector  x  is  augmented  by  3u,  quantifying  the  former,  and  introducing  as 
the  control  signal  Au,  each  defined  respectively  by  the  following  equations: 

Aw{ /,.)  =  «(*, )-»(/,_,) 

<5w  (ti+l)  =  d  «(t,.)  +  A  u{tt)  (23) 
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In  so  doing,  both  du  and  An  carry  weights  in  the  cost  equation,  constraining  both 
appropriately  with  respect  to  system  performance  requirements. 

4.3.2,  Cost  Assignment 

The  X,  U,  and  S  matrices  quantify  costs  for  state,  control,  and  cross-terms, 
respectively,  in  the  infinite  horizon  cost  equation  J,  introduced  in  Chapter  2.  As 
described  above,  the  state  vector  is  augmented  by  du,  so  X  is  a  4  by  4  matrix,  declared  in 
software  in  four  parts.  The  first  declared,  Xu,  is  a  3  by  3  matrix  and  defines  costs  for  the 
system  states.  Since  the  pull-in  voltage  was  found  by  FEM  to  be  18.36  volts,  the  cost 
assigned  is  0.0025  (the  inverse  of  the  square).  Similarly,  the  maximum  deflection  is  2 
microns,  so  the  cost  is  set  as  2.5x10* ’.  Last,  the  velocity  cost  is  set  at  unity,  with  the  goal 
of  minimizing  velocity  directly.  X12,  as  a  scalar  cross-term,  sets  the  cost  weight  for 
between-sample  dynamics  of  a  continuous  plant  not  captured  by  a  discrete  observer  and 
controller  combination.  Since  the  samples  are  taken  35  times  faster  than  the 
characteristic  time  constant  of  the  unforced  system  (modeled  as  a  damped  harmonic 
oscillator),  the  consequences  of  not  considering  interstitial  time  periods  in  the  cost 
minimization  is  assumed  to  be  negligible,  and  X12  is  concordantly  set  to  zero.  X22  is  the 
cost  weight  on  control  deviations  from  the  nominal  control  uo  and  is  a  scalar.  In  getGc() 
(described  in  the  next  section),  this  is  a  function  of  uo,  effectively  setting  X22  as  the 
inverse  square  of  20%  of  uq.  Similarly,  the  cost  weight  of  the  control  pseudorate  An  is 
quantified  by  U  and  is  defined  initially  as  lxlO4,  to  model  a  maximum  current  output  of 
the  controller  is  10  mA  (later  tuned  to  1x10  to  allow  a  more  robust  response).  Finally,  S 
is  set  to  zero  using  the  same  reasoning  as  X12. 
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4.3.3.  Dlqr  to  Get  Gc 

Having  augmented  the  state  transition  matrix  and  discrete  control  injection 
matrix  Bd  to  reflect  the  state  augmentation  described  in  the  previous  section,  and  defined 
the  cost  matrices  X,  U,  and  S,  getGc()  then  calls  the  MATLAB  script  for  discrete  linear 
quadratic  regulator  calculation,  dlqr(),  passing  each  of  these  matrices  as  inputs.  dlqr(), 
for  Discrete  Linear  Quadratic  Regulator,  calculates  regulator  gains  that  minimize  the 
quadratic  cost  function  specified  earlier,  as  well  as  the  infinite  horizon  solution  of  the 
discrete-time  Riccati  equation.  Using  dlqr()  here  implies  the  assumption  of  certainty 
equivalence;  that  is,  contributions  from  noise  sources  are  not  taken  into  account  in  the 
regulation  problem.  Assuming  the  inputs  pass  dlqr  requirements  for  a  closed  solution  (<1> 
and  Bd  must  be  stabilizable;  X  and  U  positive  definite;  and  no  unobservable  modes  in  the 
unit  circle),  getGc()  outputs  the  dlqr()  solution  for  Gc. 

4.3.4.  Dynamics  Embedding  via  II,  Kx  and 

With  regulator  gains  in  hand,  effort  is  now  spent  ensuring  that  the  controlled 
variable  yc  drives  the  equilibrium  to  the  desired  setpoint  yy 

(0-7)  Bd  =ro" 

C  Dy\  LvJ  (24) 

Assuming  the  dimensions  of  u  and  yc  are  equal,  and  the  left-hand  matrix  is  invertible,  II 
is  defined  [62]  so  as  to  solve  for  an  equilibrium  solution  as  a  function  of  yy. 

■j01J(®-i)  nAToIJn,,  n,,T o 

.«„J  L  c  nj  LvJ  Ln,  n,,J_vJ  (25) 

Perturbation  variables  may  now  be  defined  as  deviations  from  these  equilibrium  values: 
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(26) 


&  =  x(/I)-x0  =  x(?1)-ni2yd 
Su  =  »(6)  -  »0  =  “(6)-  n22^d 
4’c  =^(0-^ 

These  perturbation  variables  are  passed  to  the  cost  equation  /,  and  X,  U,  and  S  contain 
weights  appropriate  for  perturbation,  rather  than  for  the  full  state  values.  In  this  way,  the 
canonical  LQ  regulator  equation,  introduced  in  Chapter  2,  is  defined  in  terms  of 
perturbation  variables,  in  which  the  definitions  of  the  perturbation  variables  can  be 
substituted  to  produce  a  type-0,  nonzero  setpoint  controller  [62]: 

»((.)  =  »o  -Gc(r,)[x(f,)-Xo] 

=  -Gc(t,  )x(0  +  [n22  +Gc(t,)n12]yd  (27) 

As  a  type-0  controller,  this  equation  ensures  that  the  system  will  settle  to  an  equilibrium, 
though  not  necessarily  with  zero  steady-state  error,  i.e.,  to  vy;  inevitable  modeling  errors 
cause  imprecision  in  II,  creating  steady-state  errors.  As  such,  a  tenn  proportional  to  the 
total  regulation  error  is  motivated.  It  can  be  shown  that  embedding  the  system  dynamics 
into  the  controller  gains,  and  specifically  adding  a  signal  to  the  perturbation  controller 
proportional  to  the  regulation  error,  provides  the  type-1  characteristic  necessary  to 
accommodate  regulation  errors  [62]: 

K  y.]=[oa  gj 

=  1)-A'.[v((;)-.r(/;  l)]+/.'.[v1(/  l-vJ(/  ,)]  (29) 

4.3.5.  Description  of  Control  Flow 


nu  n 12 
n„  n„ 


(28) 
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A  block  diagram  of  the  controller  is  shown  in  Figure  15.  Before  the  initial  time  to, 
the  system  is  assumed  to  have  reached  steady  state;  here,  the  system  starts  at  quiescence. 
At  to,  Vd  changes  from  yd, old  to  yv/,  which  sets  the  initial  control  signal  proportional  to  the 
regulation  error,  as  the  first  and  second  terms  of  Equation  (24)  are  zero  at  to.  Using  the 
propagate  and  measurement  processes  described  above,  the  Kalman  Filter  produces  a 
refined  state  update  (and  residuals  for  diagnostics  as  motivated  in  Section  2.8),  which  is 
passed  to  the  deterministic  PI  controller.  The  controller  revises  the  control  signal  and 
applies  the  signal  simultaneously  to  a  sample-and-hold  memory  buffer  and  the  FBMD 
itself  (the  latter  only  being  modeled  via  get_meas()  and  geterr()).  At  the  start  of  the  next 
sample  period,  the  process  repeats  itself. 

4.4.  Script  Implementation 

This  section  serves  as  the  end-to-end  description  of  the  algorithm  (see  Appendix), 
tying  together  subroutines  detailed  in  previous  sections.  The  main  script  begins  by 

Disturbances 


Figure  15:  Controller  block  diagram  assuming  certainty  equivalence  [62] 
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declaring  physical  constants,  sampling  frequency,  and  sample  vector  length.  State  and 
measurement  noise  covariance  magnitudes  are  set  next,  and  are  followed  by  translation 
matrices  H  and  C,  which  map  states  to  measurement  space,  z,  and  controlled  variable 
space,  yc.  The  desired  deflection,  yd,  is  declared  next  and  used  to  calculate  Myd)  and 
u{yd).  These  latter  values  are  used  as  the  linearization  nominal  for  the  Jacobian  matrices 
F  and  B,  which  are  discretized  by  the  sub-script  equiv_discrete()  to  produce  outputs  <1>, 
Bd,  and  Q(1.  As  described  above,  the  system  dynamics  model,  in  the  form  of  <1>  and  Bd,  is 
embedded  in  controller  gains  via  calculation  of  II.  The  subroutine  getGc()  is  called  next 
and  outputs  controller  gains  associated  with  the  infinite  horizon  linear  quadratic 
regulation  problem,  as  described  above  and  as  solved  by  MATLAB’s  dlqr()  routine.  The 
final  controller  gains,  Kx  and  K=,  are  calculated  using  the  fonnulas  above. 

MATLAB  suffers  speed  penalties  for  variable-sized  arrays,  such  as  those  that 
expand  with  every  sample.  To  mitigate  this  penalty,  x,  P,  z,  u,  and  the  residual  matrices 
are  preloaded  as  zero  matrices  of  length  equal  to  the  sample  history  length,  as  declared  at 
the  beginning  of  script.  Last  declarations  are  the  filter  initial  conditions:  the  state  vector 
x  is  assumed  to  start  from  quiescence;  the  first  control  signal  is  calculated  using  the 
control  law  derived  above;  and  the  initial  covariance  matrix  P  is  found  by  trial  and  error 
to  be  10*1. 

The  Kalman  Filter  proper,  here  instantiated  as  a  for  loop,  starts  by  propagating  the 
state  vector  and  covariance  matrix  one  sample  forward.  This  is  followed  by  the  “update” 
process  begins  by  calculating  the  denominator  of  the  Kalman  gain,  A,  which  will  be  used 
later  in  the  script  for  visualization  of  the  residuals.  After  the  Kalman  gain  K,  a  scalar,  is 
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calculated,  get_meas()  is  called  to  obtain  the  discrete  measurement  z,  and  from  it  the  state 
estimate  mapped  to  measurement  space  (the  product  of  H  and  x,  effectively  retaining  the 
first  state)  is  subtracted  to  produce  the  residual.  This  residual — the  difference  of  the 
estimated  and  measured  voltage  across  the  FBMD — is  scaled  by  the  Kalman  gain  and 
added  to  the  state  estimate  to  produce  the  “updated”  state  estimate.  The  covariance 
matrix  is  updated  as  well  by  subtracting  the  product  of  K ,  H,  and  P  from  P.  The  last  step 
in  the  for  loop  is  calculating  the  deviation  of  covariance  by  taking  the  square  root  of  the 
diagonal  of  P.  This  for  loop  repeats  for  the  entire  length  of  the  time  history  declared  at 
the  beginning  of  the  script. 

Upon  completion  of  the  for  loop,  a  matrix  of  error  values  quantifying  the 
difference  between  truth,  as  defined  above,  and  the  observer’s  estimate  for  each  sample 
period  is  initialized  and  calculated  using  the  geterr()  subroutine.  Finally,  the  script 
produces  two  figures  as  visualizations  of  perfonnance.  The  first  figure  contains  two 
subplots  with  the  time  history  of  voltage  and  deflection  errors,  respectively,  and  deviation 
bounds.  The  second  figure  plots  the  time  history  of  residuals  and  confidence  bounds. 
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5.  Results  and  Conclusions 


5.1.  Introduction 

This  chapter  presents  and  discusses  the  results  of  applying  the  KF  and  LQ  controller 
developed  in  Chapter  4  to  the  FBMD  introduced  in  Chapter  1  and  analyzed  in  Chapters  2 
and  3.  The  first  section  discusses  the  KF  performance,  to  include  analysis  of  the 
residuals  and  errors  from  the  truth  model.  The  next  section  presents  the  successful 
control  of  FBMD  past  one-third  of  the  gap.  Conclusions  follow  this  section,  and  the 
chapter  finishes  with  suggestions  for  future  work. 

5.2.  Observer  Results 

The  linearized  KF  presented  in  Chapter  4  was  first  employed  with  various  constant 
control  inputs,  and  then  with  the  LQ  controller  also  developed  in  Chapter  4.  The 
modeling  software  suites  themselves  imposed  the  largest  difficulty.  COMSOL  does  not 
do  well  in  simulating  simple  circuits  in  conjunction  with  a  multiphysics  model,  and  just 
the  model  alone  did  not  allow  for  the  simulation  of  circuital  charge  flow  as  a  result  of 
conservation  of  energy.  With  no  current,  all  COMSOL  simulations,  herein  used  as  truth 
data,  were  conducted  with  the  mirror  plate  and  ground  held  at  constant  potential.  The 
key  point  is  that  the  MATLAB  script,  in  simulating  the  observer  design,  allows  for  a 
varying  potential  across  the  FBMD  with  a  constant  Norton-equivalent  applied  current, 
whereas  the  COMSOL  simulations  are  exactly  opposite.  The  end  result  is  that,  while  the 
filter  is  tuned  to  minimize  residuals  and  maintain  a  residual  mean  around  zero,  the  error 
with  respect  to  “truth”  results  are  less  definitive. 
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Filter  performance  is  tuned  primarily  through  adjustment  of  the  measurement  noise 
covariance  and  secondarily  through  adjustment  of  the  filter  noise  covariance.  The  first 
simulation  takes  as  input  a  constant  current  of  10  pA,  which  results  in  10  volts  steady- 
state  across  the  FBMD  and  a  deflection  of  97  mn.  The  observer  is  linearized  around  the 
analytical  equilibrium  state  and  control  values.  Figure  16  shows  that  the  state  estimates 
are  accurate  and  well-behaved,  while  Figure  17  demonstrates  the  residuals  are  steady, 
small  compared  state  values,  and  centered  about  zero. 


)eflection  (m) 


Figure  16:  Voltage  across  FBMD  in  blue  (V)  and  deflection  in  green  (m)  vs.  time  (s) 
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Residual  and  Uncertainty  (2  c) 


Figure  17:  Residuals  and  filter-calculated  uncertainty  in  dashed  lines  vs.  time;  constant 
input  at  10  uA 
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Figure  18:  Top  graph  shows  observer  error  in  estimate  of  voltage;  bottom  graph  error  in 
estimate  of  deflection.  Constant  input  of  10  uA.  Diverging  voltage  error  is  an  anomaly 
of  get_err(). 
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The  second  simulation  shown  takes  as  input  a  constant  current  of  1 1 .29  pA, 
which  analytically  corresponds  to  an  equilibrium  deflection  of  1.5  pm.  Figure  19 
demonstrates  the  snap-in  very  clearly,  as  the  deflection  (shown  as  the  green  line) 
increases  asymptotically.  Residuals  and  filter-computed  uncertainty  in  Figure  20  are 
increased  on  average  over  the  previous  case,  but  remain  steady  and  centered  about  zero; 
this  implies  that,  despite  being  in  an  unstable  operating  point,  the  linearized  model  holds 
as  valid.  The  error  plots  in  Figure  21  demonstrate  the  shortcoming  of  the  geterr()  script; 
since  the  COMSOL  simulation  produced  no  deflections  of  1.5  pm,  the  lack  of  data 
manifests  as  accumulating  error  and  diverging  covariance. 


Figure  19:  KF  estimate  of  voltage  across  FBMD  (blue,  in  V)  and  rapidly  diverging 
deflection  (green,  in  m)  vs.  time  (s)  for  constant  input  current  of  1 1 .29  pA 
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Residual  and  Uncertainty  (2  o) 


Figure  20:  KF-calculated  residuals  and  uncertainty  for  constant  input  of  1 1.29  pA 
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Figure  21:  Errors  in  estimate  of  voltage  across  FBMD  (top)  and  deflection  (bottom)  vs 
time.  Note  the  divergence  in  deflection  confidence,  corresponding  to  the  divergence  of 
the  system. 
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5.3.  Controller  Results 


With  a  working  linearized  KF,  the  LQ  controller  is  implemented  and  analyzed.  The 
following  set  of  figures  show  reaching  and  maintaining  desired  deflections  of  1 .0  pm  and 
1.5  pm,  respectively — each  well  beyond  one-third  of  the  gap. 

Figure  21  shows  the  deflection  and  control  signal  resulting  from  setting  the  desired 
deflection  yd  to  1  micron.  The  mirror  deflection  is  smooth  and  stable,  with  no  overshoot 


Figure  22:  Voltage  across  FBMD  (V,  in  blue)  and  deflection  (m,  in  green)  vs.  time.  The 
FBMD  achieves  a  steady  deflection  of  1  pm. 
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Residual  and  Uncertainty  (2  o) 
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Figure  23:  Residuals  and  fdter-computed  uncertainty  for  yd=1.0  pm,  tuned 
liberally  (i.e.,  Q  set  “optimistically”  or  less  than  actual  model  error). 

and  a  rise  time  of  1  millisecond.  Interestingly,  the  controller  changes  polarity  just  before 
the  deflection  reaches  the  point  of  instability  at  one-third  of  the  quiescent  gap,  go.  Figure 
22  shows  the  residual  data  for  the  first  observer  state  (voltage  across  the  FBMD)  and 
filter-computed  covariance.  As  delineated  in  Section  2.8,  model  adequacy  is  reflected  by 
the  steady  covariance,  lack  of  bias  in  the  residuals,  and  low  root-mean-square  value 
relative  to  the  magnitude  of  the  steady-state  control  signal.  However,  the  covariance 
curves  fail  to  encapsulate  the  residual  data,  but  rather  are  less  than  the  absolute  peak 
values  of  the  residual  data.  This  implies  that  the  KF  is  underestimating  the  variance  of 
the  data,  which  may  be  mitigated  by  increasing  the  dynamics  noise  strength  Q. 
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Figure  24:  Voltage  across  FBMD  (V,  in  blue)  and  deflection  (m,  in  green) 
vs.  time  (sec)  for  yd=1.5  (am 

For  the  case  in  which  the  desired  deflection  yd  is  set  to  1.5  micron,  the  resulting 
control  signal,  deflection,  residual  data,  and  covariance  are  shown  as  before  in  Figure  23 
and  Figure  24.  Figure  23  shows  successful  deflection  through  three-quarters  of  the  gap, 
smoothly  reaching  equilibrium  in  0.2  milliseconds  without  overshoot.  As  before,  Figure 
24  shows  that  the  filter  is  tuned  optimistically  and  underestimates  the  covariance.  A  key 
difference  between  Figure  24  and  Figure  22,  however,  is  that  the  residuals  start  much 
higher,  but  as  equilibrium  is  established,  the  residuals  reduce  to  a  level  closer  to  the  KF’s 
estimate. 
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Residual  and  Uncertainty  (2  o) 


Figure  25:  Residuals  and  fdter-computed  uncertainty  for  yd=1.5  pm,  tuned 
liberally. 

5.4.  Conclusions 

The  KF-based,  LQ  control  implementation  for  an  FBMD  application  is  shown  to  be 
able  to  achieve  and  maintain  deflections  beyond  the  theoretical  snap-in  point  at  one  third 
of  the  gap.  The  implementation  here  is  well  suited  for  digital  implementation,  and 
eminently  extendable  to  an  array  of  FBMDs.  Deflection  inference  via  electrical 
measurement  has  been  successfully  demonstrated  as  a  faster  and  simpler  way  to  realize 
feedback  control.  Furthermore,  precalculation  of  controller  gains  reduce  required 
computational  power  for  embedded  control;  one  might  imagine  a  small  set  of  desired 
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deflection  points  for  a  particular  set-and-hold  application,  and  associated  controller  gains 
precalculated  as  a  look-up  table. 

5.5.  Future  Work 

Future  work  may  include  a  number  of  directions.  First  and  foremost,  simulation  of 
feedback  in  a  comprehensive,  multi-physics  environment  is  required  to  fully  understand 
any  cross-coupling  between  the  FBMD  and  control  circuitry.  Second,  and  perhaps  in  lieu 
of  such  an  environment,  an  experiment  may  be  set  up  such  that  a  deflected  beam  be 
steered  to  a  distance  that  would  require  FBMD  deflection  of  more  than  one-third  of  the 
gap.  A  modem,  commercial  microcontroller  may  well  be  able  to  read  voltage  across  the 
FBMD,  source  a  small  current,  and  implement  in  real-time  the  precalculated  control  law 
presented  in  Chapter  4. 

The  LQ  control  law  may  be  extended  to  a  command  generated  tracker  (CGT)  control 
law,  with  which  the  deflection  may  follow  a  trajectory  more  complicated  than  a  simple 
set-and-hold.  With  CGT  control,  the  FBMD  may  be  capable  of  “painting,”  or  tracking,  a 
target  in  real-time,  rather  than  discrete  pointing  control.  Other  types  of  MEMS  actuation 
may  be  considered  as  well,  such  as  electrothermal  or  microfluidic.  Controller 
development  for  other  phenomena  would  proceed  in  much  the  same  way:  by  describing 
the  basic  physics,  analyzing  the  system  for  deviations,  creating  a  model,  and  translating 
that  model  into  state  space.  The  more  precise  the  model,  the  higher  performance  the  KF 
is  capable  of  achieving. 
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Appendix 


%Start  main  script 

%Physical  Constant  declaration 

Rint=le6;  %source  resistance 

e0=8 . 854187817e-12 ;  %free  space  permittivity 

a=le-8./4;  %effective  plate  area,  accounting  for  fringe  fields 
eOA=eO .  *a; 

m=4 . 66e-ll . /4 ;  %calculated  plate  mass 

ks=15 . 059 . /4 ;  %simulated  effective  spring  constant 

g0=2e-6;  %initial  gap  space 

b=1.1885e-5;  %fluid  damping  coefficient 

Vpi=sqrt ( ( 8 . *ks . *g0 . A3) . / (27 . *eOA) ) ;  %analytical  pull-in  voltage 

tc=sqrt (m. /ks) ;  %system  time  constant 

dt=50e-9;  %nanosecond  samples,  1GHz  DAq 

t_last=2500e-6; 

t=0 : dt : t_last; 

sam_last=length (t)  ; 


%Filter  Model  parameters 
ql=le-2 ; 
q2=le-8 ; 
q3=l; 

Rw=0 .06; 

R  =  (Rw) ;  %  Measurement  noise  covariance 

Q  =  [ql  0  0;0  q2  0;0  0  q3];  %  Dynamics  noise  strength 


%Mapping 
H  =  [1  0  0  ]  ; 

C  =  [0  1  0  ]  ; 

%Initialize  model  state  and  control  signal 
yd=9 .6996e-8;%le-6; 

u0= (gO-yd) . *sqrt ( (2 . *ks . *yd) .  / (eOA. *Rint . A2) ) ; 
x0= [uO . *Rint ; yd; 0 ] ; 

%  Discretize  state  dynamics  matrix  using  Jacobian 

F= [ (xO (2) -gO) . / (Rint . *eOA)  (xO ( 1 ) -uO . *Rint) . / (Rint . *eOA)  0; 

0  0  1; 

eOA. *x0 (1) . * (gO-xO (2) ) . A-2 . /m  -ks . /m+e0A. *x0 (1) .  A2  .  * (gO- 
xO  (2 ) )  .  A-3  . /m  -  b./m]; 

B=[ (gO-xO (2) ) . /eOA; 0 ; 0 ] ; 

G=eye (3) ; 

[phi,  Bd  , Qd] =equiv_discrete (F, B, G,  Q,  dt)  ;  %EENG  765, Lt  Col  Vazquez 
Gd=eye (3); 

Pi=inv ( [phi -eye (3)  ,  Bd;C,  0]  )  ; 

Gc=getGc (phi , Bd, uO ) ;  %Get  LQ  gains 

Kx=Gc ( 1 : 3) *Pi ( 1 : 3, 1 : 3) +Gc ( 4 ) *Pi ( 4 , 1 : 3) ;  %Calc  prop,  gain 
Kxi=Gc ( 1 : 3) *Pi ( 1 : 3, 4 ) +Gc (4 ) *Pi ( 4 , 4 ) ;  %Calc  pseudo-rate  gain 
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%Preloads  for  speed 
x=zeros (3, sam_last)  ; 

P=zeros (3,3, sam_last) ; 
z=zeros (2, sam_last) ; 
u=zeros ( 1 , sam_last)  ; 

A=zeros ( 1 , sam_last) ; 
sigma_f=zeros (3, sam_last) ; 
residual=zeros ( 1 , sam_last) ; 
e=zeros ( 2 , sam  last); 


%  Filter  Initial  conditions 
x  ( : , 1)  =  [ le-15  le-15  le-15]; 

u (1) =Kxi* (yd-C*x ( : , 1) ) ;  %Initial  control  signal  assuming  u0=0  and  x 
steady  state 
P ( : , : , 1) =eye (3) .*lel; 


%  Kalman  filter 
for  k=2 : sam_last 
%  Propagate 

x ( : , k) =  phi*x(:, (k-1) ) +Bd*u (k-1 ) ; 

P ( : ,  : , k) =  phi*P ( : ,  : ,  (k-1) ) *phi ' +Gd*Qd*Gd'  ; 

%  Update 

A ( : ,  k) =  H*P  (  : ,  : ,  k) *H '  +  R; 

K=P  (  : ,  :,k)*H'*(A(:,k)A-l); 

z ( : , k) =get_meas ((k-1) . *dt, k . *dt , [ z ( 1 , ( k-1 ) ) , z (2 , (k-l)),0],u(k- 
1 ) ) ; % (tO , tf , xO, u) ; 

residual  (:,k)=z(l,k)-H*x(:,k)  ; 

x ( : ,  k) =x ( : ,  k) +K* residual (:,k);  % (3x1) +  (3xl) *(lxl) 
u(k)=u(k-l)-Kx* (x ( : , k) -x ( : , k-1) ) +Kxi* (yd-C*x ( : , k-1 ) ) ; 

P ( : ,  : , k) =P ( : ,  : , k) -K*H*P ( : ,  : ,  k)  ; 
sigma_f ( : , k) =sqrt (diag ( P ( :  ,  : ,  k)  )  )  ; 
end  %  End  time  loop 


%  Compute  error  states 

e=geterr (t, x (1, : ) , x (2, : ) , u) ;  %Generate  truth  from  FEM  data 

%  Compute  statistics  of  the  error  states  and  plot 
figure ( 1 ) 
for  j  =1 : 2 

subplot ( 2 , 1 ,  j  )  ; 

plot (t, e ( j , : ) , t, 2*sigma_f ( j , : ) , ' r-- ' , t, -2*sigma_f ( j , : ) , ' r-- ' ) 

end 

f igure ( 1 ), subplot (21 1 ); grid  on 

title ([' State  Error  and  Uncertainty  (2\sigma) ' ] ) 
ylabel ( ' Voltage ' )  ; 

subplot ( 212 ) ; ylabel ( ' Gap ' ) ; grid  on 
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figure (2 ) 

plot (t, residual, ' o- ' , t, 2 . *sqrt (A),'r--',t,-2. *sqrt (A) , ' r-- ' ) 

title ([' Residual  and  Uncertainty  (2\sigma)']) 

end 


function  [Gc]  =  getGc (phi , Bd, uO ) 
phia=[phi  Bd;0  0  0  1]; 

Bda= [zeros (3,1) ; 1] ; 

%should  be  incremental ! 

Xll= [0.0025  0  0 ; 0  2.5ell  0;0  0  l];%3e3];  %max  V=20  g=2e-6  v=0.0185 
X12=zeros (3,1); 

X22  =  l . /  (0 .2 . *u0)  . A2; 

Xcost= [Xll  X 1 2 ; X 1 2 '  X2 2 ] ; 

Ucost=le3 ; %le4 ;  %assume  source  max  10mA 
S=zeros (4,1) ; 

%S= [1; 0; 0]  ; 

%E=eye(2);  %same  dimension  as  Bd 
% [Kc, L, Gc, report ]  =  dare (phi, Bd, Xcost, Ucost) 

[Gc, S, E] =dlqr (phia, Bda, Xcost, Ucost, S) 
end 

function  out  =  get_meas (tO, tf , xO, u) 

%calculates  analytical  solution  of  state  dynamics  given  start  time 
to, 

%stop  time  tf,  previously  analytically  calc'd  state  as  initial 
condition 

%x0,  and  control  input  u 
%Output  is  [voltage,  gap] 

[t, y] =ode45 (Spotent, [tO  tf ] , xO,  [  ] , u) ; 

out=y (size (y, 1) , 1 : 2) ' ;  %just  take  the  last  voltage  and  gap  solutions 

%R_nois=0 . 01 . *randn;  %emulate  +/-10%  resistor  tolerance 

%g_nois=0 . 1 . *randn;  %optical  sensor  error 

R_nois=0 . 1 . *  rand- . 05 ; 

g_nois=0 . 1 . *  rand- . 05 ; 

out ( 1 , 1 ) =out (1,1) . * ( l+R_nois) ;  %voltage 
out ( 2 , 1 ) =out (2,1) . * ( l+g_nois) ;  %gap 
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function  dxdt=potent (t, x, u) 

Rint=le6 ; %source  resistance 

e0=8.854187817e-12; 

a=le-8 . /4 ; 

eOA=eO .  *a; 

m=4 . 66e-ll . /4; 

ks=l 5 . 05  9 . /4 ; 

g0=2e-6; 

b=9.4779e-4; 

dxdt  ( 3 )  =  0 . 5  .  *eOA .  *x  ( 1 )  .  A2  .  *  (gO-x  (2  )  )  . A-2  .  /  m-b  .  *x  ( 3 )  .  /m-ks  .  *x  ( 2  )  .  /m; 
dxdt ( 2 ) =x ( 3 ) ; 

dxdt ( 1 ) = ( (gO-x (2 ) ) . /eOA) . * (u-x ( 1 ) . /Rint ) ; 
dxdt=dxdt ' ; 
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function  err=geterr (tobs, Vobs, Dobs, u) 

load  leap. mat 

load  Qfem_CxV.mat 

load  dispfem.mat 

stop=length (tobs) ; 

Qobs  =  zeros ( 1 , stop)  ; 
err=zeros (2, stop) ; 

Qobs  ( 1 ) =0 ; 
for  i=2 : stop 

Qobs (i) =trapz (tobs(l,l:i), (u-Vobs (1,1:1) . /Rint ) ) ; 
for  j =2 : length (Qfem)  %find  1st  time  Qfem>Qobs 
[r, c] =f ind ( (Qf em ( j ) . q-Qobs (i) ) >0, 1) ; 

if (isempty (r) )%find( (Qfem ( j ) . q-Qobs (i))>0,l, 'first'))) 
elseif(Qfem(j) . q==Qobs ( i) )  Qf lag=l 

else  % index ( 1, i) =f ind ( (Qfem(i) . q-Qobs) <0,1, 'last') 
index=j ; 

%output=Qf em ( j ) . v ( 1 ) ; 
break 

end 

end 

%err (1, i) =interpl (Qfem ( index) . q, Qfem ( index) . v, Qobs ,' spline ' ) -Vobs ; 
%now  linearly  interpolate  between  voltages 

if (r==l  &&  index==2)  %when  Qobs  is  less  than  all  Qfem  data 
last=length (Qfem ( index- 1 ) . v) ; 

m= (Qfem ( index) .q(r)-0) ./ (Qfem (index) . v (r ) -0) ; 

err ( 1 , i ) = ( (Qobs (i) -Qfem ( index) . q ( r ) +m . *Qf em ( index) . v ( r ) ) . /m) - 
Vobs  (i) ; 

elseif (r==l  &&  index>2)  %when  Qobs  is  btn  discrete  FEM  steps 
last=length (Qfem ( index- 1 ) . v) ; 

m= (Qfem ( index) . q ( r ) -Qfem ( index- 1) . q ( last ) ) ./ (Qfem ( index) . v ( r ) - 
Qfem (index-1) .v (last) ) ; 

err ( 1 , i ) = ( (Qobs (i) -Qfem ( index) . q ( r ) +m . *Qf em ( index) . v ( r ) ) . /m) - 
Vobs  (i) ; 

else  %when  Qobs  falls  within  an  FEM  sim 
last=r-l ; 

m= (Qfem ( index) . q ( r ) -Qfem ( index) . q ( last ) ) . / (Qfem ( index) . v ( r ) - 
Qfem (index) . v (last) ) ; 

err ( 1 , i ) = ( (Qobs (i) -Qfem ( index) . q ( r ) +m . *Qf em ( index) . v ( r ) ) . /m) - 
Vobs  (i) ; 
end 
end 

%err=Qobs ; 

for  i=l : length (Dispfem) 

disp_resam (i ) . ts=re sample (Dispfem (i ) . ts, tobs) ; 

end 

for  tau=l : length (tobs)  %for  each  trajectory  time  tau 
for  j=l:51  %and  each  input 

gvsv (tau) . d ( j , 2 ) =disp_resam ( j ) . ts . data (tau) ; 

gvsv (tau) . d ( j , 1 ) =Qf em ( j ) . v ( 1) ; 

end 

%tau; 


73 


err(2,tau)=interpl (gvsv (tau) . d ( : , 1 ) , gvsv (tau) . d ( : , 2 ) , Vobs ( tau) , ' spli 
ne ' ) -Dobs (tau) ; 

%errd (tau) =temp (tau) -x (2, tau) ; 

end 

err(2,tau)=interpl (gvsv (tau) . d ( : , 1 ) , gvsv (tau) . d ( : , 2 ) , Vobs ( tau) , metho 
d, ' spline ' ) ; 
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